¿Cómo evaluaron las computadoras Apollo funciones trascendentales como seno, arcotangente, logaritmo?

Dado que el código del Apolo 11 está en GitHub, pude encontrar el código que parece una implementación de funciones de seno y coseno: vea aquí el módulo de comando y aquí el módulo de aterrizaje lunar (parece que es el mismo código).

Para mayor comodidad, aquí hay una copia del código:

 # 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.

El comentario

# SINGLE PRECISION SINE AND COSINE

indica que lo siguiente es de hecho una implementación de las funciones seno y coseno.

La información sobre el tipo de ensamblador utilizado se puede encontrar en Wikipedia .

Explicación parcial del código:

La subrutina SPSINrealmente calcula$\sin(\pi x)$, y SPCOScalcula$\cos(\pi x)$.

La subrutina SPCOSprimero suma un medio a la entrada y luego procede a calcular el seno (esto es válido debido a$\cos(\pi x) = \sin(\pi (x+\tfrac12))$). El argumento se duplica al comienzo de la SPTsubrutina. Es por eso que ahora tenemos que calcular$\sin(\tfrac\pi2 y)$por$y=2x$.

La subrutina POLLEYcalcula una aproximación polinomial casi de Taylor de$\sin(\tfrac\pi2 x)$. Primero, almacenamos$x^2$en el registro SQ (donde$x$denota la entrada). Esto se utiliza para calcular el polinomio.

que se parecen a los primeros coeficientes de Taylor para la función$\frac12 \sin(\tfrac\pi2 x)$.

¡Estos valores no son exactos! Así que esta es una aproximación polinomial, que está muy cerca de la aproximación de Taylor, pero aún mejor (ver más abajo, también gracias a @uhoh y @zch).

Finalmente, el resultado se duplica con el DDOUBLcomando y la subrutina POLLEYdevuelve una aproximación a$\sin(\tfrac\pi2 x)$.

En cuanto a la escala (primero la mitad, luego el doble, ...), @Christopher mencionó en los comentarios que el número de punto fijo de 16 bits solo podía almacenar valores de -1 a +1. Por lo tanto, es necesario escalar. Consulte aquí para obtener una fuente y más detalles sobre la representación de datos. Los detalles de las instrucciones del ensamblador se pueden encontrar en la misma página.

¿Qué tan precisa es esta aproximación casi de Taylor? Aquí puede ver un gráfico en WolframAlpha para el seno, y parece una buena aproximación para$x$desde$-0.6$a$+.6$. La función coseno y su aproximación se representan aquí . (Espero que nunca hayan tenido que calcular el coseno para un valor$\geq \tfrac\pi2$, porque entonces el error sería desagradablemente grande).

@uhoh escribió un código de Python , que compara los coeficientes$C_{1/2}, C_{3/2}, C_{5/2}$del código Apollo con los coeficientes de Taylor y calcula los coeficientes óptimos (basados ​​en el error máximo para$-\tfrac\pi2 \leq x \leq \tfrac\pi2$y error cuadrático en ese dominio). Muestra que los coeficientes de Apolo están más cerca de los coeficientes óptimos que los coeficientes de Taylor.

En este gráfico las diferencias entre$\sin(\pi x)$y se muestran las aproximaciones (Apollo/Taylor). Se puede ver que la aproximación de Taylor es mucho peor para$x\geq .3$, pero mucho mejor para$x\leq .1$. Matemáticamente, esto no es una gran sorpresa, porque las aproximaciones de Taylor solo se definen localmente y, por lo tanto, a menudo solo son útiles cerca de un solo punto (aquí$x=0$).

Tenga en cuenta que para esta aproximación polinomial solo necesita cuatro multiplicaciones y dos sumas ( MPy ADen el código). Para la computadora de orientación Apollo , la memoria y los ciclos de CPU solo estaban disponibles en pequeñas cantidades.

Hay algunas formas de aumentar la precisión y el rango de entrada, que habrían estado disponibles para ellos, pero daría como resultado más código y más tiempo de cálculo. Por ejemplo, explotar la simetría y la periodicidad del seno y el coseno, usar la expansión de Taylor para el coseno o simplemente agregar más términos de la expansión de Taylor habría mejorado la precisión y también habría permitido valores de entrada grandes arbitrarios.

También pediste el logaritmo, así que hagamos esto también. A diferencia de las funciones seno y coseno, esta no se implementa únicamente con un enfoque similar a una serie de Taylor. El algoritmo se basa en cambiar la entrada y contar el número de cambios necesarios para llegar a la escala requerida. No sé el nombre de este algoritmo, esta respuesta en SO describe el principio básico.

La LOGimplementación es parte del módulo CGM_GEOMETRY y está etiquetada

SUBRUTINA PARA CALCULAR EL LOGOTIPO NATURAL

La rutina utiliza la NORMinstrucción ensambladora, que, según su documentación , desplaza a la izquierda el número en el acumulador (registro "MPAC"), hasta que termina con un valor$\geq {0.5}$y "casi$1$" ^[1] , y escribe el número de operaciones de desplazamiento realizadas en la ubicación de memoria especificada como su segundo argumento (el significado matemático de la operación de desplazamiento a la izquierda es exponenciación binaria$2^n$, los exponentes en el argumento de un logaritmo se pueden expresar como factores y los productos en el argumento se pueden expresar como sumas, por lo que la simplificación de$ln(2^c \cdot scaled)$en$c \cdot const + ln(scaled)$trabaja, donde$c$es el número de turnos y$const$es el valor precalculado de$ln(2)$).

Luego usa un polinomio de tercer grado para aproximar$ln$en ese intervalo, con coeficientes codificados ^[3] :

Eventualmente, el conteo de turnos obtenido anteriormente se vuelve a multiplicar (por la constante 0.0216608494 ^[2] , usando SHORTMP).

La presión de optimización debe haber sido tan alta que no corrigieron el signo invertido antes de regresar de la subrutina, lo que requirió que todos los sitios de llamadas lo tuvieran en cuenta.

Aplicación de la subrutina logarítmica:
por ejemplo, como parte de la predicción de rango en el control de reentrada.

^[1] el formato de almacenamiento para un número de doble precisión se construyó a partir de dos palabras de 16 bits donde el MSB de cada uno es el signo y forma una representación del rango en complemento a uno$-1 < x < 1$pero el LSB es un bit de paridad. Por lo tanto, tratamos con un formato de 30 bits que contiene dos bits de signo, lo que genera un dolor de cabeza para la implementación del emulador.

^[2] el lenguaje ensamblador ACG lo permite CLOG2/32como nombre identificador. Esto causó más dolores de cabeza para la implementación del emulador .

^[3] ¿cómo se encontraron los coeficientes? Los comentarios del código sobre la lista de ensamblaje del módulo de trigonometría interpretativa (sí, los astronautas podrían hacer que el ACG interprete instrucciones dinámicas) sugieren que el método se basó en trabajos de C. Hastings, especialmente Aproximaciones para computadoras digitales. Prensa de la Universidad de Princeton, 1955 . El polinomio más complejo de ese tipo en ACG es uno de séptimo orden, mismo módulo, para calcular$arccos(x) \approx \sqrt{1-x} \cdot P(x)$).

Una adición a la respuesta de @supinf :

a) La inicial DOUBLEenSPT

b) desbordamientos para la entrada x (registro A) por encima de +0,5 (+90°) y subdesbordamientos para x por debajo de -0,5 (-90°). En este caso, Aes >+1 o <-1 y lo siguiente TSalmacena el valor corregido (efectivamente suma uno, si está por debajo de -1, o resta uno, si está por encima de +1) TEMKy lo establece Aen +1 ( desbordamiento) o -1 (desbordamiento). Además , se ignora el salto TCFa .POLLEY

c) se XCHintercambia TEMKcon A, por lo que Acontiene el valor corregido y TEMK±1 ahora.

d) INDEXsuma el valor de TEMK(±1) al valor de la siguiente ADinstrucción, que corrige silenciosamente los desbordamientos. Como LIMITSes igual a -0,5, esto da como resultado una suma de 0,5 en el caso de desbordamiento (-0,5 + +1 = 0,5) y una resta de 0,5 en el caso de subdesbordamiento (-0,5 + -1 = -1,5 = -0,5) .

e) COMniega el valor de A– esto incluye invertir el bit de desbordamiento – y

f) el final ADsuma uno en el caso de overflow y resta uno en el caso de underflow. ADno realiza una corrección de desbordamiento antes de la adición y establece el indicador de desbordamiento después. Entonces, cada valor desbordado (>+135° y <-135°) volverá al rango [-1,+1].

Si esto se ADdesborda o se desborda (no veo forma de que esto suceda), se establece Aen ±1, salta ARG90y se establece Aen -( LIMITS+ A), que es -(-0.5+1)=-(+0.5)=-0.5 en el caso de desbordamiento y -(-0.5-1)=-(-1.5)=-(-0.5)=+0.5 en el caso de subdesbordamiento. Inicialmente pensé que esto sucedería para x > +135° o x < -135°, pero no parece ser el caso.

Pero esto de ajustar la carcasa <-90° y >+90° me parece un poco incorrecto. Yo esperaría que la línea f fuera de (+0.5,+1.0) a (+1.0,+0.0) y de (-0.5,-1.0) a (-1.0,-0.0). Ese sería el caso si COMsigue directamente XCHsin el paso d.

Corríjame, si tengo algunas piezas incorrectas, recientemente vi ese código y traté de resolverlo con el conjunto de instrucciones AGC .