Wie haben die Apollo-Computer transzendente Funktionen wie Sinus, Arkustangens, Logarithmus ausgewertet?

Da sich der Apollo 11-Code auf GitHub befindet, konnte ich den Code finden, der wie eine Implementierung von Sinus- und Cosinus-Funktionen aussieht: Siehe hier für das Befehlsmodul und hier für den Mondlander (es sieht so aus, als wäre es derselbe Code).

Der Einfachheit halber hier eine Kopie des Codes:

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

Der Kommentar

# SINGLE PRECISION SINE AND COSINE

zeigt an, dass das Folgende tatsächlich eine Implementierung der Sinus- und Kosinusfunktionen ist.

Informationen über den verwendeten Assemblertyp finden Sie auf Wikipedia .

Teilerklärung des Codes:

Das Unterprogramm SPSINrechnet tatsächlich$\sin(\pi x)$, und SPCOSrechnet$\cos(\pi x)$.

Die Subroutine SPCOSaddiert zuerst eine Hälfte zur Eingabe und fährt dann mit der Berechnung des Sinus fort (dies ist gültig wegen$\cos(\pi x) = \sin(\pi (x+\tfrac12))$). Das Argument wird am Anfang des SPTUnterprogramms verdoppelt. Deshalb müssen wir jetzt rechnen$\sin(\tfrac\pi2 y)$Pro$y=2x$.

Die Subroutine POLLEYberechnet eine nahezu taylorische polynomische Approximation von$\sin(\tfrac\pi2 x)$. Zuerst speichern wir$x^2$im Register SQ (wobei$x$bezeichnet den Eingang). Daraus wird das Polynom berechnet

die wie die ersten Taylor-Koeffizienten für die Funktion aussehen$\frac12 \sin(\tfrac\pi2 x)$.

Diese Werte sind nicht exakt! Das ist also eine Polynom-Näherung, die der Taylor-Näherung sehr nahe kommt, aber noch besser ist (siehe unten, auch dank @uhoh und @zch).

Schließlich wird das Ergebnis mit dem DDOUBLBefehl verdoppelt, und die Subroutine POLLEYgibt eine Annäherung an zurück$\sin(\tfrac\pi2 x)$.

Bezüglich der Skalierung (erst halbieren, dann verdoppeln, ...) erwähnte @Christopher in den Kommentaren, dass die 16-Bit-Festkommazahl nur Werte von -1 bis +1 speichern könne. Daher ist eine Skalierung erforderlich. Hier finden Sie eine Quelle und weitere Details zur Datendarstellung. Details zur Montageanleitung finden Sie auf der gleichen Seite.

Wie genau ist diese Fast-Taylor-Näherung? Hier sehen Sie ein Diagramm auf WolframAlpha für den Sinus, und es sieht nach einer guten Annäherung für aus$x$von$-0.6$zu$+.6$. Hier ist die Kosinusfunktion und ihre Näherung aufgetragen . (Ich hoffe, sie mussten nie den Kosinus für einen Wert berechnen$\geq \tfrac\pi2$, denn dann wäre der Fehler unangenehm groß.)

@uhoh hat einen Python-Code geschrieben , der die Koeffizienten vergleicht$C_{1/2}, C_{3/2}, C_{5/2}$aus dem Apollo-Code mit den Taylor-Koeffizienten und berechnet die optimalen Koeffizienten (basierend auf dem maximalen Fehler für$-\tfrac\pi2 \leq x \leq \tfrac\pi2$und quadratischer Fehler in diesem Bereich). Es zeigt, dass die Apollo-Koeffizienten näher an den optimalen Koeffizienten liegen als die Taylor-Koeffizienten.

In diesem Diagramm die Unterschiede zwischen$\sin(\pi x)$und die Annäherungen (Apollo/Taylor) werden angezeigt. Man sieht, dass die Taylor-Näherung viel schlechter ist$x\geq .3$, aber viel besser für$x\leq .1$. Mathematisch ist dies keine große Überraschung, da Taylor-Approximationen nur lokal definiert sind und daher oft nur in der Nähe eines einzelnen Punktes (hier$x=0$).

Beachten Sie, dass Sie für diese Polynomnäherung nur vier Multiplikationen und zwei Additionen benötigen ( MPund ADim Code). Für den Apollo Guidance Computer standen Speicher und CPU-Zyklen nur in geringer Stückzahl zur Verfügung.

Es gibt einige Möglichkeiten, die Genauigkeit und den Eingabebereich zu erhöhen, die für sie verfügbar gewesen wären, aber dies würde zu mehr Code und mehr Rechenzeit führen. Beispielsweise hätte das Ausnutzen der Symmetrie und Periodizität von Sinus und Cosinus, die Verwendung der Taylor-Entwicklung für Cosinus oder das einfache Hinzufügen weiterer Terme der Taylor-Entwicklung die Genauigkeit verbessert und auch beliebig große Eingabewerte ermöglicht.

Sie haben auch nach dem Logarithmus gefragt, also machen wir das auch. Im Gegensatz zu den Sinus- und Kosinusfunktionen wird diese nicht nur mit einem Taylor-Reihen-ähnlichen Ansatz implementiert. Der Algorithmus basiert auf dem Verschieben der Eingabe und dem Zählen der Anzahl der Verschiebungen, die erforderlich sind, um den erforderlichen Maßstab zu erreichen. Ich kenne den Namen dieses Algorithmus nicht, diese Antwort auf SO beschreibt das Grundprinzip.

Die LOGImplementierung ist Teil des CGM_GEOMETRY -Moduls und gekennzeichnet

SUBROUTINE ZUR BERECHNUNG DES NATÜRLICHEN LOGS

Die Routine verwendet die NORMAssembler-Anweisung, die laut ihrer Dokumentation die Zahl im Akkumulator ("MPAC"-Register) nach links verschiebt, bis sie mit einem Wert endet$\geq {0.5}$und „fast$1$" ^[1] , und schreibt die Anzahl der durchgeführten Schiebeoperationen in die als zweites Argument angegebene Speicherstelle (die mathematische Bedeutung der linken Schiebeoperation ist binäre Potenzierung$2^n$, Exponenten im Argument eines Logarithmus können als Faktoren ausgedrückt werden und Produkte im Argument können als Additionen ausgedrückt werden, also die Vereinfachung von$ln(2^c \cdot scaled)$hinein$c \cdot const + ln(scaled)$funktioniert, wo$c$ist die Schichtzahl und$const$ist der vorberechnete Wert von$ln(2)$).

Dann verwendet es ein Polynom dritten Grades zur Annäherung$ln$in diesem Intervall mit fest codierten Koeffizienten ^[3] :

Schließlich wird die zuvor erhaltene Verschiebungszahl wieder multipliziert (mal die Konstante 0,0216608494 ^[2] , unter Verwendung von SHORTMP).

Der Optimierungsdruck muss so hoch gewesen sein, dass sie das umgekehrte Vorzeichen nicht behoben haben, bevor sie aus der Subroutine zurückkehrten, sondern von allen Aufrufstellen verlangten, dies stattdessen zu berücksichtigen.

Anwendung des Logarithmus-Unterprogramms:
zum Beispiel als Teil der Reichweitenvorhersage in der Wiedereintrittssteuerung.

^[1] Das Speicherformat für eine Zahl mit doppelter Genauigkeit wurde aus zwei 16-Bit-Wörtern aufgebaut, wobei das MSB von jedem das Vorzeichen ist und eine Einerkomplementdarstellung des Bereichs bildet$-1 < x < 1$aber das LSB ist ein Paritätsbit. Wir haben es also mit einem 30-Bit-Format zu tun, das zwei Vorzeichenbits enthält, was einige Kopfschmerzen bei der Emulatorimplementierung verursacht.

^[2] die ACG-Assemblersprache lässt CLOG2/32als Identifikatornamen zu. Dies verursachte einige weitere Kopfschmerzen bei der Implementierung des Emulators .

^[3] Wie wurden die Koeffizienten gefunden? Code- Kommentare zur Montageliste des interpretativen Trignonometriemoduls (ja, Astronauten konnten das ACG dazu bringen, dynamische Anweisungen zu interpretieren) legen nahe, dass die Methode auf Arbeiten von C. Hastings basierte, insbesondere Approximations for Digital Computers. Princeton University Press, 1955 . Das komplexeste Polynom dieser Art in ACG ist eines der siebten Ordnung, gleicher Modul, der zu berechnen ist$arccos(x) \approx \sqrt{1-x} \cdot P(x)$).

Eine Ergänzung zur Antwort von @supinf :

a) Die Initiale DOUBLEinSPT

b) Überläufe für Eingang x (Register A) über +0,5 (+90°) und Unterläufe für x unter -0,5 (-90°). In diesem Fall Aist >+1 oder <-1 und im Folgenden wird TSder korrigierte Wert gespeichert (effektiv eins addiert, wenn er unter -1 liegt, oder eins subtrahiert, wenn er über +1 liegt) TEMKund Aauf +1 gesetzt ( Überlauf) oder -1 (Unterlauf). Außerdem wird der Sprung TCFzu POLLEYignoriert.

c) XCHvertauscht TEMKmit A, Aenthält also jetzt den korrigierten Wert und TEMK±1.

d) INDEXAddiert den Wert von TEMK(±1) zum Wert der nächsten ADAnweisung, die stillschweigend Überläufe korrigiert. Da LIMITSgleich -0,5 ist, ergibt sich im Überlauffall eine Addition von 0,5 (-0,5 + +1 = 0,5) und im Unterlauffall eine Subtraktion von 0,5 (-0,5 + -1 = -1,5 = -0,5) .

e) COMnegiert den Wert von A– dazu gehört auch das Invertieren des Overflow-Bits – und

f) das Finale ADaddiert eins im Überlauffall und subtrahiert eins im Unterlauffall. ADführt vor der Addition keine Überlaufkorrektur durch und setzt danach das Überlauf-Flag. Jeder übergelaufene Wert (>+135° und <-135°) kommt also wieder in den Bereich [-1,+1].

Wenn dies ADunter/überläuft (ich sehe keine Möglichkeit, dass dies passieren könnte), wird es Aauf ±1 gesetzt, springt auf ARG90und wird Aauf -( LIMITS+ A) gesetzt, was -(-0,5+1)=-(+0,5)=-0,5 Zoll ist im Überlauffall und -(-0,5-1)=-(-1,5)=-(-0,5)=+0,5 im Unterlauffall. Ich dachte zunächst, dies würde für x > +135° oder x < -135° passieren, aber das scheint nicht der Fall zu sein.

Aber diese Anpassung des Gehäuses <-90° und>+90° sieht für mich irgendwie falsch aus. Ich würde erwarten, dass die Linie f von (+0,5, +1,0) bis (+1,0, +0,0) und von (-0,5, -1,0) bis (-1,0, -0,0) reicht. Das wäre der Fall, wenn COMdirekt XCHohne Schritt d folgt.

Bitte korrigieren Sie mich, wenn ich einige Teile falsch mache, ich habe diesen Code erst kürzlich gesehen und versucht, ihn mit dem AGC-Befehlssatz herauszufinden .