162

Navigation with a sextant or maneuvers using gimbal angles might be two examples of cases where an Apollo computer might need to do trigonometry.

Trigonometric functions like sine, arctangent, etc. are transcendental, as are logarithms. You can't evaluate these functions with a simple expression built on multiplication or division for example, without at least an iterative algorithm.

An engineer on the ground would grab a slide rule for two or three digits, and for more go to a book of trig, log and other tables for more digits. Between two lines you could interpolate by hand for even more digits.

But how did the Apollo computers evaluate transcendental functions, or how at least were calculations that required the use of transcendental functions implemented in the programs?

Glorfindel
  • 723
  • 1
  • 9
  • 22
uhoh
  • 148,791
  • 53
  • 476
  • 1,473
  • 7
    Are you sure about the need to do trigonometry for gimbaling? I'm pretty sure it would be easier to stay within vector math if you use a linear actuator. Anyway trignonometry functions are usually implemented by table lookup or approximated as a taylor expansion. No source directly related to Apollo, sorry. – Christoph Sep 27 '18 at 11:14
  • 1
    ...or both, the taylor expansion preparing lookup tables on startup, in case permanent storage is less abundant than RAM. – SF. Sep 27 '18 at 11:16
  • @Christoph First sentence starts with two instruments that likely produce angular data. Also, what was done in the 1960's in computers in space, is not covered by how things are done now. I'll add the history tag to make that even clearer. You might also consider where tables would have to be stored, there was precious little memory in the Apollo computers. – uhoh Sep 27 '18 at 11:16
  • @Christoph have a look here: https://youtu.be/9YA7X5we8ng?t=1283 and also here https://youtu.be/YIBhPsyYCiM?t=273 – uhoh Sep 27 '18 at 11:30
  • 14
    @Christoph: Already in 1956 we had a better algorithm than Taylor, namely CORDIC. I don't know if that was used in Apollo. – MSalters Sep 27 '18 at 11:43
  • 1
    Storage (permanent and dynamic both) was at a huge premium so it's not surprising they didn't keep a table. – Russell Borogove Sep 27 '18 at 16:30
  • 1
    One could ask this same question of e.g. calculators, or really any computer. – NeutronStar Sep 27 '18 at 19:44
  • 3
    @Joshua consider that these were developed in the early 1960's. as well as considering the size and weight constraints for landing on the moon and returning to Earth, this isn't just "any computer". In fact, technologies developed for computers in the space program helped pave the way for "personal" scientific calculators a decade later. It's the totality of the situation that makes this particular question compelling. – uhoh Sep 27 '18 at 20:19
  • 7
    @MSalters: CORDIC needs tables of constants, about 50 values. Not very usefull for the Apollo computers with core rope memory for program and constants. – Uwe Sep 27 '18 at 21:03
  • @uhoh, I agree the totality of the situation is compelling, but the meat of the answer (Taylor expansions) is the same as if I asked how my computer does it, at least how I understand these calculations happen today. – NeutronStar Sep 28 '18 at 03:12
  • 1
    @Joshua while modern implementations may include one or more taylor expansions as a seed (depending on which function, sin is easy, found here), that's just the beginning of how modern computers do double precision transcendentals. What I learned from this answer is exactly how it was done in this case, the degree of resulting precision (~1E-04) they had to work with, the effort that went into pre- and post-scaling and why, and the spartan coding. – uhoh Sep 28 '18 at 05:13
  • @Joshua all of these are unique to the specific set of constraints on this particular, one(few)-of-a-kind computer. To your reductionist use of "is the same as" the only thing I can say is "no it isn't", to which you could reply "yes it is" and we could continue ad infinitum – uhoh Sep 28 '18 at 05:19
  • 7
    @Joshua exactly - every digital computer does it by this or similar methods. Another way to look at it is "where did the tables in the books we used to look this stuff up in come from?" In the old days, some poor schmucks had to grind through those series by hand, and later with the aid of mechanical calculators that could multiply and divide. The general field of study is called "Numerical methods" and you can easily find numerous books under that subject heading. – Jamie Hanrahan Sep 28 '18 at 06:51
  • @JamieHanrahan there's SE site too! Coincidentally I've just asked What algorithm does (did) Excel use for Bessel functions that is discontinuous at x=8? – uhoh Sep 28 '18 at 06:54
  • 1
    @Uwe actually CORDIC was developed specifically for digital navigation computers in aircraft in 1956. – Chris Stratton Sep 28 '18 at 17:45
  • for log it's even possible to output the value to an analog computer, compute it and read the result back in an ADC. Those things can be done very quickly in an electric analog computer. I don't know if it works trigonometry but a mechanical analog computer can do that albeit slower – phuclv Sep 29 '18 at 07:50
  • 1
    @Cris Straton: From wikipedia: "when a hardware multiplier is available (e.g., in a DSP microprocessor), table-lookup methods and power series are generally faster than CORDIC". The Apollo computer had a pretty fast multiplication, Memory Cycle time: 11.7 microseconds. Addition Time: 23.4 microseconds. Multiplication Time: 46.8 microseconds. – Uwe Sep 29 '18 at 12:04
  • 2
    @phuclv: there were analog electronic circuits of that time for logarithm, exponentiation, adding, multiplication and square root but not for trigonometric functions like sin, cos and tan. – Uwe Sep 29 '18 at 13:57

3 Answers3

302

Since the Apollo 11 code is on GitHub, I was able to find the code that looks like an implementation of sine and cosine functions: see here for the command module and here for the lunar lander (it looks like it is the same code).

For convenience, here is a copy of the code:

 # Page 1102
            BLOCK   02

# SINGLE PRECISION SINE AND COSINE

            COUNT*  $$/INTER
SPCOS       AD      HALF            # ARGUMENTS SCALED AT PI
SPSIN       TS      TEMK
            TCF     SPT
            CS      TEMK
SPT         DOUBLE
            TS      TEMK
            TCF     POLLEY
            XCH     TEMK
            INDEX   TEMK
            AD      LIMITS
            COM
            AD      TEMK
            TS      TEMK
            TCF     POLLEY
            TCF     ARG90
POLLEY      EXTEND
            MP      TEMK
            TS      SQ
            EXTEND
            MP      C5/2
            AD      C3/2
            EXTEND
            MP      SQ
            AD      C1/2
            EXTEND
            MP      TEMK
            DDOUBL
            TS      TEMK
            TC      Q
ARG90       INDEX   A
            CS      LIMITS
            TC      Q       # RESULT SCALED AT 1.

The comment

# SINGLE PRECISION SINE AND COSINE

indicates, that the following is indeed an implementation of the sine and cosine functions.

Information about the type of assembler used, can be found on Wikipedia.

Partial explanation of the code:

The subroutine SPSIN actually calculates $\sin(\pi x)$, and SPCOS calculates $\cos(\pi x)$.

The subroutine SPCOS first adds one half to the input, and then proceeds to calculate the sine (this is valid because of $\cos(\pi x) = \sin(\pi (x+\tfrac12))$). The argument is doubled at the beginning of the SPT subroutine. That is why we now have to calculate $\sin(\tfrac\pi2 y)$ for $y=2x$.

The subroutine POLLEY calculates an almost Taylor polynomial approximation of $\sin(\tfrac\pi2 x)$. First, we store $x^2$ in the register SQ (where $x$ denotes the input). This is used to calculate the polynomial $$ ((( C_{5/2} x^2 ) + C_{3/2} ) x^2 + C_{1/2}) x. $$ The values for the constants can be found in the same GitHub repository and are

$$\begin{aligned} C_{5/2} &= .0363551 \approx \left(\frac\pi2\right)^5 \cdot \frac1{2\cdot 5!}\\ C_{3/2} &= -.3216147 \approx -\left(\frac\pi2\right)^3 \cdot \frac1{2\cdot 3!}\\ C_{1/2} &= .7853134 \approx \frac\pi2 \cdot \frac12\\ \end{aligned}$$

which look like the first Taylor coefficients for the function $\frac12 \sin(\tfrac\pi2 x)$.

These values are not exact! So this is a polynomial approximation, which is very close to the Taylor approximation, but even better (see below, also thanks to @uhoh and @zch).

Finally, the result is doubled with the DDOUBL command, and the subroutine POLLEY returns an approximation to $\sin(\tfrac\pi2 x)$.

As for the scaling (first halve, then double, ...), @Christopher mentioned in the comments, that the 16-bit fixed-point number could only store values from -1 to +1. Therefore, scaling is necessary. See here for a source and further details on the data representation. Details for the assembler instructions can be found on the same page.

How accurate is this almost-Taylor approximation? Here you can see a plot on WolframAlpha for the sine, and it looks like a good approximation for $x$ from $-0.6$ to $+.6$. The cosine function and its approximation is plotted here. (I hope they never had to calculate the cosine for a value $\geq \tfrac\pi2$, because then the error would be unpleasantly large.)

@uhoh wrote some Python code, which compares the coefficients $C_{1/2}, C_{3/2}, C_{5/2}$ from the Apollo code with the Taylor coefficients and calculates the optimal coefficients (based on the maximal error for $-\tfrac\pi2 \leq x \leq \tfrac\pi2$ and quadratic error on that domain). It shows that the Apollo coefficients are closer to the optimal coefficients than the Taylor coefficients.

In this plot the differences between $\sin(\pi x)$ and the approximations (Apollo/Taylor) is displayed. One can see, that the Taylor approximation is much worse for $x\geq .3$, but much better for $x\leq .1$. Mathematically, this is not a huge surprise, because Taylor approximations are only locally defined, and therefore they are often only useful close to a single point (here $x=0$).

Note that for this polynomial approximation you only need four multiplications and two additions (MP and AD in the code). For the Apollo Guidance Computer, memory and CPU cycles were only available in small numbers.

There are some ways to increase accuracy and input range, which would have been available for them, but it would result in more code and more computation time. For example, exploiting symmetry and periodicity of sine and cosine, using the Taylor expansion for cosine, or simply adding more terms of the Taylor expansion would have improved the accuracy and would also have allowed for arbitrary large input values.

Nathan Tuggy
  • 4,566
  • 5
  • 34
  • 44
supinf
  • 2,223
  • 1
  • 8
  • 7
  • 2
    I am reading into it, but it is not easy because i am not familiar with the programming language. At first glance it looks like a taylor approximation, but i am not sure. I will edit when i have more. – supinf Sep 27 '18 at 11:33
  • 1
    I'm on it... i think it has to do with scaling. – supinf Sep 27 '18 at 12:35
  • 4
    The scaling is due to the fact that single precision could only store value from -1 to +1 :). Precision is around 13 bits which fits the single precision type (16 bits one of which is the sign bit and one was a parity bit not accessable to software). – Christoph Sep 27 '18 at 12:45
  • @Christoph So I assume that it was a fixed point representation and not floating point, correct? Also, do you have a link that i can add to the answer (or you can edit)? – supinf Sep 27 '18 at 13:02
  • 1
    @uhoh: The scaling issues are fixed now, and i included plots to compare the function with the approximation. It seems to fit for some values. – supinf Sep 27 '18 at 13:07
  • 1
    I found that information in the Virtual AGC Programmer's Manual. – Christoph Sep 27 '18 at 13:13
  • 13
    On topic: error graph. Maximum error within the domain would be around 0.0001 – SF. Sep 27 '18 at 13:13
  • 20
    "I hope they never had to calculate the cosine for a value ≥ π/2." It is not necessary using the relation cos(π - α) = -cos(α). Using this and similar relations, only the range 0 ≤ α ≤ π/4 has to be computed. These transformations may be used for sin, cos, tan and cot. There may be other functions within the Apollo 11 code using these transformations when bigger arguments are possible. – Uwe Sep 27 '18 at 15:55
  • 2
    Might some of the instructions near the start have been truncating 2x to the range +/-1 and then inverting if that wraps? If so, the behavior of cos() would have been clean for any angle. – supercat Sep 27 '18 at 18:19
  • 43
    @MagicOctopusUrn The Brogrammers among us who claim that "girls can't code" should take note that this code was developed by a woman. – Philipp Sep 28 '18 at 11:43
  • 4
    Taylor is based on derivatives in a single point, it is not a best (minimal max-error) fit over a range. So it is expected that similar, but different polynomial would be used. – zch Sep 28 '18 at 13:23
  • 1
    @supinf indeed it seems the coefficients are optimized to give good results over ±π/2. A quick check gives nearly the same numbers as in the program, rather than the nominal values of a proper taylor expansion at x=0. https://pastebin.com/UnVudQs4 Very insightful answer by the way, thank you again! – uhoh Sep 29 '18 at 11:54
  • 7
    @Philip: This code was developed by a team of mostly men and few women led by Magaret Hamilton. She did not tell here if she was the only women in the team. – Uwe Sep 29 '18 at 14:21
  • 1
    @uhoh So what they did was much better than Taylor approximation (@zch suggested this already). I have edited and also linked to your python code. Great Analysis - it shows that mean square error was more important to them than the maximal error. – supinf Sep 30 '18 at 03:37
  • 1
    @supinf the pastebin (rhymes with waste bin) script was only a quick test. Minimization can be surprisingly tricky and I didn't take the time to look into it thoroughly, so I've walked your discussion about the results back somewhat. Thanks! – uhoh Sep 30 '18 at 04:34
  • 2
    The choice of $C$ as a constant is a moderate hint that the polynomial in question isn't a Taylor expansion, but rather a quadrature expansion in Chebyshev polynomials or something similar. – E.P. Oct 01 '18 at 14:16
  • @E.P.: I hadn't thought of that, but I'd been curious how the function seemed to have such a wide useful range, since my own experience with Taylor series suggested that the useful range is much smaller. Accepting a loss of precision at points nearer zero, however, allows points further out to be much more accurate. – supercat Oct 01 '18 at 18:01
  • @uhoh Since the minimization problems are convex and there is also the message es'Optimization terminated successfully.' for resmrs and resmax, i was fairly confident in the results. Also maxiter and maxfev where not reached. But at least it is clear that the Apollo coefficients are much better than the Taylor coefficients – supinf Oct 03 '18 at 21:30
  • @supinf algorithms are not perfect, and we can have 100% faith in their perfection until we are bitten once The optimize module offers a large number of very different algorithms to choose from, I think the next step would be to try several and see just how close their results are to each other, or in this case to test it analytically to see if it's at least a local maximum. Finite size of $x$ array could make it non-convex in a quirky way, the minimized function could be non-smooth (discontinuous gradient). – uhoh Oct 04 '18 at 03:01
  • 1
    @Uwe This page credits authors of individual AGC subprograms - while male-leaning, I see 5 likely-female names on the list. https://www.ibiblio.org/apollo/ – Russell Borogove Oct 09 '18 at 03:51
  • Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on [meta], or in [chat]. Comments continuing discussion may be removed. – called2voyage Mar 20 '23 at 18:21
41

You also asked for the logarithm, so let's do this as well. As opposed to the sine and cosine functions, this one is not implemented with a Taylor series-like approach only. The algorithm is based on shifting the input and counting the number of shifts needed to arrive at the required scale. I don't know the name of this algorithm, this answer on SO describes the basic principle.

The LOG implementation is part of the CGM_GEOMETRY module, and labeled

SUBROUTINE TO COMPUTE THE NATURAL LOG

The routine uses the NORM assembly instruction, which, according to its documentation left-shifts the number in the accumulator ("MPAC" register), until it ends up with a value $\geq {0.5}$ and "nearly $1$"[1], and writes the number of shift operations performed into the memory location specified as its second argument (the mathematical meaning of the left shift operation is binary exponentiation $2^n$, exponents in the argument of a logarithm can be expressed as factors and products in the argument can be expressed as additions, so the simplification of $ln(2^c \cdot scaled)$ into $c \cdot const + ln(scaled)$ works, where $c$ is the shift count and $const$ is the precomputed value of $ln(2)$).

Then it uses a polynomial of third degree to approximate $ln$ in that interval, with hardcoded coefficients[3]:

enter image description here

Eventually the shift count obtained earlier is multiplied back in (times the 0.0216608494 constant[2], using SHORTMP).

Optimization pressure must have been so high that they did not fix the inverted sign before returning from the subroutine, requiring all call sites to take that into account instead.

Application of the logarithm subroutine:
for example as a part of range prediction in re-entry control.

---

[1] the storage format for a double precision number was built from two 16 bit words where the MSB of each is the sign and forms a one's complement representation of the range $-1 < x < 1$ but the LSB is a parity bit. So we deal with a 30 bit format containing two sign bits, causing some headache for emulator implementation.

[2] the ACG assembly language permits CLOG2/32 as an identifier name. This caused some more headache for emulator implementation.

[3] how were the coefficients found? Code comments on the assembly listing of the interpretive trignonometry module (yes, astronauts could make the ACG interpret dynamic instructions) suggest that the method was based on works by C. Hastings, especially Approximations for Digital Computers. Princeton University Press, 1955. The most complex polynomial of that kind in ACG is one of seventh order, same module, to compute $arccos(x) \approx \sqrt{1-x} \cdot P(x)$).

No Nonsense
  • 1
  • 2
  • 8
  • 21
dlatikay
  • 713
  • 4
  • 6
  • 2
    I suspect the equation ln(arg) = ln((2^a)b) = aln(2) + ln(b) is used. a is the shift count needed to get 0.5 < b < 1 so that arg = (2^a)b. But there should be a constant ln(2) = 0.693147 which I could not find. But ln(2)/32 = 0.02166 is close to the constant .0215738898. 32 .0215738898 is 0,69036. – Uwe Oct 01 '18 at 14:03
  • 2
    on second sight, there is a polynomial involved. line 274 evaluates a polynomial of the third degree (line 276, 2+1) with coefficients .031335467 etc. that follow. I'll emulate that and add a plot if I have time later. – dlatikay Oct 01 '18 at 14:24
  • 3
    @Uwe The constant for $\ln(2)/32$ is present in line 302. – Litho Oct 03 '18 at 08:27
  • 3
    0.0216608494 * 32 is 0.693147181 indeed. I'm still struggling with that weird double-precision one's complement-with-overflow and parity storage, don't get anywhere near yet with an attempted emulation in C. – dlatikay Oct 03 '18 at 10:33
  • 4
    @uhoh It looks to me that $0.031335467(1-x)+0.0130145859(1-x)^2+0.0215738898(1-x)^3$ is a good approximation of $-\ln(x)/32$ for $x$ between $0.5$ and $1$. – Litho Oct 03 '18 at 15:11
  • @Litho: good find in line 302. I should have searched more at the bottom. It is amazing that all 9 digits of the constant are the same as calcaluated with my old hp32S. A 30 bit constant may hold up to 9 decimal digits, so only the last decimal digit 4 is lost in the binary representation. – Uwe Oct 03 '18 at 20:30
  • @Litho yep, nice! – uhoh Oct 04 '18 at 03:37
  • @dlatikay I've got another question about AGC code if you're interested. https://space.stackexchange.com/q/31196/195 – Russell Borogove Oct 09 '18 at 03:52
3

One addition to the answer by @supinf:

a) The initial DOUBLE in SPT

b) overflows for input x (register A) above +0.5 (+90°) and underflows for x below -0.5 (-90°). In this case, A is >+1 or <-1 and the following TS stores the corrected value (effectively add one, if it is below -1, or subtract one, if it is above +1) in TEMK and sets A to the +1 (overflow) or -1 (underflow). Furthermore the jump TCF to POLLEY is ignored.

c) XCH swaps TEMK with A, so A contains the corrected value and TEMK ±1 now.

d) INDEX adds the value from TEMK (±1) to the value of the next AD instruction, which silently corrects overflows. Because LIMITS is equal to -0.5, this results in an addition of 0.5 in the overflow case (-0.5 + +1 = 0.5) and in a subtraction of 0.5 in the underflow case (-0.5 + -1 = -1.5 = -0.5).

e) COM negates the value of A – this includes inverting the overflow bit – and

f) the final AD adds one in the overflow case and subtracts one in the underflow case. AD does not perform a overflow correction prior to addition and sets the overflow flag after. So every overflowed value (>+135° and <-135°) will come back into the range [-1,+1].

Visualisation of the calculations a-f

If this AD under/overflows (I see no way this could happen), it sets A to ±1, jumps to ARG90 and sets A to -(LIMITS+A), which is -(-0.5+1)=-(+0.5)=-0.5 in the overflow case and -(-0.5-1)=-(-1.5)=-(-0.5)=+0.5 in the underflow case. I initially thought, this would happen for x > +135° or x < -135°, but this does not seem to be the case.

But this adjusting the case <-90° and >+90° looks kinda wrong to me. I would expect the line f to be from (+0.5,+1.0) to (+1.0,+0.0) and from (-0.5,-1.0) to (-1.0,-0.0). That would be the case if COM follows directly XCH without step d.

Please correct me, if I have some pieces wrong, I just recently saw that code and tried to figure it out with the AGC Instruction Set.

uhoh
  • 148,791
  • 53
  • 476
  • 1,473
sivizius
  • 131
  • 2