Margaret Heafield Hamilton was the director of the Software Enginee...
This code was written in AGC assembly language to be run in the Apo...
The subroutine $SPCOS$ calculates cos(𝜋𝑥). It starts by adding 1/2...
The $POLLEY$ subroutine calculates a polynomial approximation of $\... The 16-bit fixed-point number used in AGC could only store values f... Apollo 11 implementation of trigonometric functions Margaret H. Hamilton March 1969 This code was submitted by Margaret H. Hamilton 1 (Programming Leader Apollo Guidance and Navigation) in March 1969 and was used to evaluate trigonometric functions essential for navigation on both the Command Mod- ule and Lunar Lander that ﬁrst landed humans on the Moon. 1 Code is available at https://github.com/chrislgarry/Apollo-11/blob/ 27e2acf88a6345e2b1064c8b006a154363937050/Luminary099/SINGLE_PRECISION_ SUBROUTINES.agc @Luís I think you might want to exchange C12 and C52. Thanks, this is the interactive graph https://www.desmos.com/calculator/fmnu5rs6d6 The 16-bit fixed-point number used in AGC could only store values from -1 to +1. Fixed, Thanks Mark! See also: http://jeanmariechauvet.com/papers/curttrig.html the factor 12 in the sin argument should be 1/2, I think Margaret Heafield Hamilton was the director of the Software Engineering Division of the MIT Instrumentation Laboratory, which developed on-board flight software for NASA's Apollo program. ![](https://upload.wikimedia.org/wikipedia/commons/thumb/a/aa/Margaret_Hamilton_in_action.jpg/1280px-Margaret_Hamilton_in_action.jpg) $0.0363551 . x^5 - 0.3216147 . x^3 + 0.7853134 . x$ Indeed the factors where switched. Using the formulae above produces the correct graph. Thank you Roger K. for the nice words. Hasting *does* in fact explain the iterative procedure he used in Chapter 3 of the book. Apparently he did use a starting set of points to fit and then refined this trial reference iteratively, in what looks like a variant of the Remez algorithm -- although I find the detailed explanations wanting. Based on this revisiting of the procedure, I've tried to answer your questions and remarks in the revised version of the paper (see specifically Addenda). Indeed,$SPSIN$does start by calling Transfer to Storage. But the point is not to move the value to memory. The point is actually to test if the value has overflowed. If it has, then the negative of the value is calculated by the "CS TEMPK" instruction. If the routine was entered from$SPCOS$and we have just added 1/2, then we could have an overflow. True sin(pi + theta) == sin(-theta), (so the code does no harm) but why does the code allow for an already overflowed value to be passed into the$SPSIN$routine? Or, put another way, why isn't$SPT$the entry point for the sine routine? I'm guessing there must be callers that rely on having overflow logic in the sine function. The$POLLEY$subroutine calculates a polynomial approximation of$\sin(\frac{\pi}{2}y)$. $$\frac{1}{2}\sin\left(\frac{\pi}{2}y\right)=((C_{5/2}y^2+C_{3/2})y^2 + C_{1/2})y$$ $$C_{1/2}= 0.0363551 \\ C_{3/2}= -0.3216147 \\ C_{5/2}= 0.7853134$$ Note that the instructions$AD$and$MP$are used for addition and multiplication and$SQ$is the register where$y^2$is stored. The terms of this approximation are close to the terms of the Taylor approximation. $$\frac{\pi^5}{7680} \approx 0.0398463 \\ -\frac{\pi^3}{96} \approx -0.3229820 \\ \frac{\pi}{4} \approx 0.7853981$$ 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$-\pi/2 < y < \pi/2$. Apparently, most of the calculations necessary during the mission were supposed to be within this region. In this graph we can see the differences between$\sin(\pi y)$and the Taylor and Apollo approximations. ![](https://i.imgur.com/4MCppAF.png) The final result is doubled, so$POLLEY$actually returns$\sin\left(\frac{\pi}{2}y\right)$. Here's a graph comparing the Sine function with the Apollo approximation. ![](https://i.imgur.com/UkWZn7K.png) You may wonder why the "/2" in the polynomial coefficients. The symbols names are C5/2, C3/2 and C1/2. C1/2 is the coefficient for the linear term; C3/2 for the cubic term; and C5/2 for the quintic term. The "/2" is part of the name and does not cause a division. Indeed, the coefficients have already been divided by 2. The reason that the coefficients have been divided by two (compared to an error minimized polynamial for sin(x * pi/2)) is that the Apollo computer number system is a 1 bit sign with 14 data bits with the implied binary point before the data bits. That is, the first data bit, when set to 1, represent 1/2, the next 1/4, etc. This number representation can represent values in the OPEN interval -1 to 1. That is, not including 1 and not including -1. The linear coefficient, before the halving, would have been 1.5706268. Since this is outside of the representable range, the code uses half of (all) the "normal" coefficients for the calculation, and then, at the end doubles the result to compensate. The other coefficients are in range for the AGC as are their halves. Luis mentioned that the Apollo team chose these coefficients rather than those from a straight, but truncated, Taylor series to minimize the error in the range. This is true, but it's not clear that the Apollo team did the choosing (or minimizing of error). Cecil Hasting, Jr., of the RAND corporation, in his 1955 book, "Approximations for Digital Computers", gives the exact same (well if you double them!) coefficients for a fifth degree polynomial for sin(x * pi/2), for x in [-1, 1]. ![Image of sin(x * pi/2) from Hastings Approximations for Digital Computers](https://i.imgur.com/SjmNrxv.jpg) The subroutine$SPCOS$calculates cos(𝜋𝑥). It starts by adding 1/2 to the input (AD operation) and then proceeds to calculate the sine, since cos(𝜋𝑥)=sin(𝜋(𝑥+1/2)). Note that all the arguments for the routines are in multiples of$\pi $. The subroutine$SPSIN$calculates the Sine function. It starts by calling the Transfer to Storage (TS) instruction, which saves the contents of the accumulator in the memory location referenced by the address field. It then proceeds to call the$SPT$routine, using the TCF (Transfer Control to Fixed) instruction. The$SPT$routine doubles the argument, which means that we are now calculating$\sin(\frac{\pi}{2}y)$for$y=2x$, and calls$POLLEY\$. This code was written in AGC assembly language to be run in the Apollo Guidance Computer (AGC). The AGC was produced for the Apollo program and was installed on board each Apollo command module (CM) and Apollo Lunar Module (LM). The AGC provided computation and electronic interfaces for guidance, navigation, and control of the spacecraft. ![](https://upload.wikimedia.org/wikipedia/commons/7/79/Agc_view.jpg) The user interface to the AGC was the DSKY (display and keyboard and usually pronounced "DIS-kee"). The DSKY had an array of indicator lights, numeric displays, and a calculator-style keyboard. Commands were entered numerically, as two-digit numbers: Verb, and Noun. Verb described the type of action to be performed and Noun specified which data were affected by the action specified by the Verb command. ![](https://upload.wikimedia.org/wikipedia/commons/thumb/f/f1/Apollo_DSKY_interface.svg/1920px-Apollo_DSKY_interface.svg.png) For a detailed overview of the AGC architecture and Operation I highly recommend reading the book: [The Apollo Guidance Computer: Architecture and Operation by Springer](https://www.amazon.com/Apollo-Guidance-Computer-Architecture-Operation/dp/1441908765). I graphed the equation of estimation of sine but it looked like this https://www.desmos.com/calculator/k5ex58z4vs what is the problem? Very nice Jean-Marie! The only think I'd add is that I believe the Hastings polynomial is the solution to minimizing the maximum RELATIVE error in the domain. And I think that makes sense if you think about how you'd use the trig functions. C.T. Fike in "Computer Evaluation of Mathematical Functions", 1968 on page 68 mentions a Chebyshev theorem on the uniqueness of polynomial approximations that amazingly (to me!) applies to both absolute and relative errors (and more!) Also, I don't think Hastings necessarily used a discrete set of points to fit, like the trial reference, but rather iteratively found the points where the relative error was largest, sought to minimize that error, and then found the new place where the relative error was the largest and so on.