I'm a fan of Margaret Hamilton, but this is fake. This code is not by Margaret Hamilton and it's not from 1969.
Margaret Hamilton has stated that her first assignment was the abort code FORGETIT in the Apollo 5 program Sunburst, and this exact single-precision sine/cosine routine already existed in Aurora [an earlier version of the software] earlier than March 1966 -- FORGETIT wasn't added to Sunburst until later that year (maybe October-November).
Source: I talked to people who have researched the Apollo code in detail. For more information on the Apollo software releases, see: https://www.ibiblio.org/apollo/Luminary.html
I don’t read a claim that she wrote that code, only that she submitted it somewhere.
I read that as some indication as to the correctness of the claim that this is what was used (a bit like the difference between an email from Linus claiming something about early Linux and a thread on Reddit where X claims his neighbor knows somebody whose nephew worked with Linus, and had this on an old floppy, believing it to apply to early Linux)
Everything I have read about Margaret Hamilton makes me happy; such a wonderful engineer/developer, so seeing her name on this as well (not unexpected though) just enforces that.
Indeed. This was a completely new architecture, and as such, there was probably very little "prior art" to go on. Maybe there were similar functions written for other machines, but this one would pretty much guarantee dead astronauts if it weren't implemented properly.
Note at the beginning it handles overflow conditions.
Funny how these pioneers wrote bulletproof raw machine code and assembly, while today we have a slew of code analysis tools and "safe" programming languages, but everything still crashes six ways from Sunday.
>The reason why the Apollo team chose these terms instead of the Taylor ones was because these were actually optimized for minimum error in the interval −π/2<y<π/2
How would these coefficients be derived in the first place?
We did this in my quantum mechanics class, 8.05x on edX. Consider the functions f_n(x) = x^n to be basis vectors in a vector space. Define the dot product of two vectors (aka functions) to be the integral of their product between -pi/2 and +pi/2. You'll see that your basis vectors aren't orthogonal. So, use Graham Schmidt to make them orthogonal, and normalize them while you're at it.
Now you just take the dot product of these new basis vectors with your target function.
Worth noting that this minimizes expected error. IMO, it's much more common to want to minimize maximum (worst-case) error, which requires an algorithm like Remez exchange. You can modify this algorithm to work with relative error (as suggested in the article) instead of absolute error.
If one understands linear algebra (and has Haskell output of matrices with nice eigenvalues), then reviewing this spelling is the only preparation needed to teach a semester course.
Also, something doesn't add up here because they were calculating sin(pi*y/2). I suspect it's optimized for -1 < y < 1, so the argument to sin is in [-pi/2, pi/2].
In a similar vein, I can't go without linking to the Fast Inverse Square Root function. Similar kind of thing, a "machine-friendly" way to approximate a different mathematical function.
Fast inverse square root is an ultra-fancy hack for a very particular processor. Apollo's trigonometric functions are the most natural optimal approximation that could be done.
Fast inverse square root isn't that ultra fancy. It just happens that the IEEE representation lines up the exponent to the left, allowing you to do the obvious multiply-the-exponent thing, but in integer mode, followed by a standard Newton-Raphson round. The only surprising part was that this was faster than calling the CPU's native floating point instructions.
I mean that modern cpus have already a fast instruction to compute that. The apollo code does not seem to be ultra-optimized, but it is a natural implementation of the polynomial approximation.
AGC assembler uses unusual mnemonics, so it can be difficult to understand for someone who doesn't know the architecture.
The CPU has onle accumulator which instructions act on. The argument is normally an address. The relevant instructions for this code snippet is:
AD - Add the value pointed to by the argument to the accumulator
TS - Transfer to storage. Write the accumulator to the argument
TCF - Transfer control to fixed memory. Jumps to an instruction in fixed storage
DOUBLE - Assembler macro compiling to AD A. I.e. double the accumulator.
XCH - This exchanges the value in the accumulator with the argument
INDEX - Adjust the address used in the next instruction by the argument
COM - Complements the accumulator. Note that the AGC uses one's complement arithmetic, so this negates the value.
EXTEND - Indicate that the next instruction is an extended instruction
MP - Multiply the accumulator with the argument
The final instruction is TS Q which jumps to the location stored in the Q register. This is the return address which is updated by the TS instruction. In other words, TS can be used both as a regular jump as well as being used as a function call, depending on what the destination does with the Q register.
I hope this helps you analyse the code further. I haven't spent much time looking at the code itself, so I am not going to attempt an analysis right now.
Single precision is a 15-bit fixed point number in the range between -1 and 1. Apparently in most cases this was enough for moon flights. There was also double and triple precision that were used sometimes.
Of course, in many places in the code you need values outside of this range and the actual number represented is scaled by some factor. However, the CPU has no internal representation of the scale factor (i.e. the exponent bits in a floating point number) so the programmer had to manually manage this.
Fermat's Library is an annotation platform, so there are annotations on the linked version - click on the circles (or rectangles, depending on whether the side panel is open) on either side of the document to see the annotations.
Each one of those bound 'books' is a complete full printout. That photo shows a bunch of different revisions / programs stacked up.
Here's one complete listing in a book printout being discussed by one of the original authors Don Eyles.
https://www.youtube.com/watch?v=H0ITFbDuJz0
If folks are interested in their own custom approximations, I can highly recommend Sollya [1].
Sollya has a great Remez implementation as well as a proven supremum norm (maximum error). I combined it with more traditional “loop over all 32-bit float values to also test my range reduction”, but being able to choose your error and number of coefficients is nice!
There is a (non-free) book on the architecture of the guidance computer. Goes into great details on the instruction set, memory layout, multitasking, virtual machine.
"The Apollo Guidance Computer: Architecture and Operation"
Margaret Hamilton has stated that her first assignment was the abort code FORGETIT in the Apollo 5 program Sunburst, and this exact single-precision sine/cosine routine already existed in Aurora [an earlier version of the software] earlier than March 1966 -- FORGETIT wasn't added to Sunburst until later that year (maybe October-November).
Source: I talked to people who have researched the Apollo code in detail. For more information on the Apollo software releases, see: https://www.ibiblio.org/apollo/Luminary.html