Как компьютеры «Аполлон» оценивали трансцендентные функции, такие как синус, арктангенс, логарифм?

Поскольку код Аполлона-11 находится на GitHub, мне удалось найти код, который выглядит как реализация функций синуса и косинуса: здесь командный модуль и здесь лунный посадочный модуль (похоже, это один и тот же код).

Для удобства привожу копию кода:

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

Комментарий

# SINGLE PRECISION SINE AND COSINE

указывает, что следующее действительно является реализацией функций синуса и косинуса.

Информацию о типе используемого ассемблера можно найти в Википедии .

Частичное объяснение кода:

Подпрограмма SPSINфактически вычисляет$\sin(\pi x)$, и SPCOSвычисляет$\cos(\pi x)$.

Подпрограмма SPCOSсначала добавляет половину ко входу, а затем переходит к вычислению синуса (это верно из-за$\cos(\pi x) = \sin(\pi (x+\tfrac12))$). Аргумент удваивается в начале SPTподпрограммы. Поэтому теперь нам нужно вычислить$\sin(\tfrac\pi2 y)$за$y=2x$.

Подпрограмма POLLEYвычисляет почти полиномиальную аппроксимацию Тейлора$\sin(\tfrac\pi2 x)$. Во-первых, мы храним$x^2$в регистре SQ (где$x$обозначает вход). Это используется для вычисления многочлена

которые выглядят как первые коэффициенты Тейлора для функции$\frac12 \sin(\tfrac\pi2 x)$.

Эти значения не являются точными! Так что это полиномиальное приближение, которое очень близко к приближению Тейлора, но даже лучше (см. ниже, также благодаря @uhoh и @zch).

Наконец, команда удваивает результат DDOUBL, и подпрограмма POLLEYвозвращает приближение к$\sin(\tfrac\pi2 x)$.

Что касается масштабирования (сначала вдвое, затем вдвое, ...), @Christopher упомянул в комментариях, что 16-битное число с фиксированной точкой может хранить только значения от -1 до +1. Поэтому необходимо масштабирование. См. здесь источник и дополнительные сведения о представлении данных. Подробности инструкций по ассемблеру можно найти на той же странице.

Насколько точно это почти Тейлоровское приближение? Здесь вы можете увидеть график на WolframAlpha для синуса, и он выглядит как хорошее приближение для$x$от$-0.6$к$+.6$. Здесь показана функция косинуса и ее аппроксимация . (Надеюсь, им никогда не приходилось вычислять косинус для значения$\geq \tfrac\pi2$, потому что тогда ошибка была бы неприятно большой.)

@uhoh написал код на Python , который сравнивает коэффициенты$C_{1/2}, C_{3/2}, C_{5/2}$из кода Аполлона с коэффициентами Тейлора и вычисляет оптимальные коэффициенты (исходя из максимальной ошибки для$-\tfrac\pi2 \leq x \leq \tfrac\pi2$и квадратичная ошибка в этой области). Это показывает, что коэффициенты Аполлона ближе к оптимальным коэффициентам, чем коэффициенты Тейлора.

На этом графике различия между$\sin(\pi x)$и отображаются приближения (Аполлон/Тейлор). Видно, что приближение Тейлора гораздо хуже для$x\geq .3$, но гораздо лучше для$x\leq .1$. С математической точки зрения это не является большой неожиданностью, потому что тейлоровские приближения определяются только локально и, следовательно, часто полезны только вблизи одной точки (здесь$x=0$).

Обратите внимание, что для этой полиномиальной аппроксимации вам нужно всего четыре умножения и два сложения ( MPи ADв коде). Для управляющего компьютера Apollo циклы памяти и процессора были доступны только в небольшом количестве.

Есть несколько способов повысить точность и диапазон ввода, которые были бы доступны для них, но это привело бы к большему количеству кода и большему времени вычислений. Например, использование симметрии и периодичности синуса и косинуса, использование разложения Тейлора для косинуса или просто добавление дополнительных членов разложения Тейлора повысило бы точность, а также позволило бы использовать произвольно большие входные значения.

Вы также просили логарифм, так что давайте сделаем и это. В отличие от функций синуса и косинуса, эта не реализуется только с помощью подхода, подобного ряду Тейлора. Алгоритм основан на сдвиге входных данных и подсчете количества сдвигов, необходимых для достижения требуемого масштаба. Я не знаю названия этого алгоритма, этот ответ на SO описывает основной принцип.

Реализация LOGявляется частью модуля CGM_GEOMETRY и помечена как

ПОДПРОГРАММА ДЛЯ ВЫЧИСЛЕНИЯ НАТУРАЛЬНОГО ЛОГА

Подпрограмма использует NORMассемблерную инструкцию, которая, согласно документации , сдвигает число в аккумуляторе (регистре "MPAC") влево до тех пор, пока не будет получено значение$\geq {0.5}$и "почти$1$" ^[1] , и записывает количество выполненных операций сдвига в ячейку памяти, указанную в качестве второго аргумента (математический смысл операции сдвига влево - двоичное возведение в степень$2^n$, показатели степени в аргументе логарифма могут быть выражены как множители, а произведения в аргументе могут быть выражены в виде сложений, поэтому упрощение$ln(2^c \cdot scaled)$в$c \cdot const + ln(scaled)$работает, где$c$это количество смен и$const$предварительно вычисленное значение$ln(2)$).

Затем он использует полином третьей степени для аппроксимации$ln$в этом интервале с жестко заданными коэффициентами ^[3] :

В конце концов полученный ранее счетчик сдвигов снова умножается (умножается на константу 0,0216608494 ^[2] , используя SHORTMP).

Давление оптимизации, должно быть, было настолько велико, что они не исправили перевернутый знак перед возвратом из подпрограммы, вместо этого требуя, чтобы все сайты вызова учитывали это.

Применение подпрограммы логарифмирования:
например, как часть предсказания дальности при повторном входе в атмосферу.

^[1] формат хранения для числа двойной точности был построен из двух 16-битных слов, где MSB каждого является знаком и формирует представление диапазона с дополнением до единицы$-1 < x < 1$но LSB - это бит четности. Таким образом, мы имеем дело с 30-битным форматом, содержащим два знаковых бита, что вызывает некоторую головную боль при реализации эмулятора.

^[2] язык ассемблера ACG позволяет CLOG2/32использовать имя идентификатора. Это вызвало дополнительную головную боль для реализации эмулятора .

^[3] как были найдены коэффициенты? Комментарии к коду ассемблерного листа интерпретирующего тригнонометрического модуля (да, астронавты могли заставить ACG интерпретировать динамические инструкции) позволяют предположить, что метод был основан на работах Ч. Гастингса, особенно Approximations for Digital Computers. Издательство Принстонского университета, 1955 г. Самый сложный многочлен такого рода в ACG — это полином седьмого порядка того же модуля для вычисления$arccos(x) \approx \sqrt{1-x} \cdot P(x)$).

Одно дополнение к ответу @supinf :

а) Начальный DOUBLEвSPT

б) переполнения для входа x (регистр A) выше +0,5 (+90°) и потери значимости для x ниже -0,5 (-90°). В этом случае Aэто >+1 или <-1, а следующее TSсохраняет исправленное значение (фактически добавьте единицу, если оно ниже -1, или вычтите единицу, если оно больше +1) TEMKи задает Aзначение +1 ( переполнение) или -1 (недополнение). TCFКроме того , переход POLLEYигнорируется.

c) XCHзаменяется TEMKна A, поэтому теперь Aсодержит исправленное значение и TEMK±1.

г) INDEXдобавляет значение из TEMK(±1) к значению следующей ADинструкции, что молча исправляет переполнения. Поскольку LIMITSравно -0,5, это приводит к прибавлению 0,5 в случае переполнения (-0,5 + +1 = 0,5) и к вычитанию 0,5 в случае недополнения (-0,5 + -1 = -1,5 = -0,5). .

e) инвертирует COMзначение A- включая инвертирование бита переполнения - и

f) финал ADдобавляет единицу в случае переполнения и вычитает единицу в случае потери значимости. ADне выполняет коррекцию переполнения перед добавлением и устанавливает флаг переполнения после. Таким образом, каждое переполненное значение (>+135° и <-135°) вернется в диапазон [-1,+1].

Если это ADпереполняется / переполняется (я не вижу, как это могло бы произойти), оно устанавливается Aв ± 1, переходит к ARG90и устанавливается Aв -( LIMITS+ A), что составляет -(-0,5+1)=-(+0,5)=-0,5 в в случае переполнения и -(-0,5-1)=-(-1,5)=-(-0,5)=+0,5 в случае недостаточного переполнения. Сначала я думал, что это произойдет при x > +135° или x < -135°, но, похоже, это не так.

Но такая настройка корпуса <-90° и >+90° мне кажется неправильной. Я ожидаю, что строка f будет от (+0,5,+1,0) до (+1,0,+0,0) и от (-0,5,-1,0) до (-1,0,-0,0). Это было бы так, если бы COMследовало непосредственно XCHбез шага d.

Пожалуйста, поправьте меня, если я ошибаюсь в некоторых частях, я только недавно видел этот код и пытался понять его с помощью AGC Instruction Set .