Изучаем FMUL из Altair Basic

nabbla1 — 16.02.2025 Со сложением с плавающей точкой чуть разобрались, как на 8080 под Altair Basic, так и под Intel 8051 (раз, два). И там, и там штука весьма и весьма мудрёная.

Посмотрим теперь, а как реализовано умножение с плавающей точкой на процессоре, который умеет только складывать и сдвигать. Как ни странно, код здесь гораздо короче, хотя исполняется, конечно, подольше (не считая особенности сложения Altair Basic, которое может затянуться на десятки миллисекунд).

Этот код определённо "с характером" - товарищ любит играться с Tail Call, разбивать код на множество отдельных функций, а тут ещё и "прерывание" применил, и впридачу "наложение команд" и самомодифицирующийся код.


Дизассемблированный код берём отсюда: https://github.com/option8/Altair-BASIC/blob/master/BASIC%20disassembly-source.lst

"Основная часть" FMUL занимает адреса от 0x08E3 до 0x0922 (всего 63 байта), но, как увидим, часть кода сидит в других местах. Начнём:

08E3: C1        FMul    POP B   ;Get lhs in BCDE
08E4: D1                POP D   ;
08E5: EF                RST 05  ; FTestSign     ;If rhs==0 then exit
08E6: C8                RZ      ;

Напоминаем: правый операнд (right-hand side, или rhs) находится в "аккумуляторе" в памяти, там зарезервировано 4 байта под это дело. А левый операнд (left-hand side, или lhs) мы передаём в стеке, причём на самой верхушке, "над" адресом возврата. Как именно это делается - пока не смотрел.

В 8080 push/pop 16-битный, поэтому первой строкой выпихиваются экспонента в регистр B, знак и старшие биты мантиссы в C, и младшие 2 байта мантиссы - в D,E.

Далее происходит то самое "прерывание", RST (Restart). Специфическая команда 8080, которая, по сути, вызывает функцию, адрес которой в данном случае 5 * 8 = 40, или 0x28. Т.е всё начало памяти "поделено" на отдельные 8-байтные блоки, и с помощью RST можно прыгать в 0x00, 0x08, 0x10, 0x18, и так далее до 0x38 = 56.

Судя по всему, в этом бейсике данные адреса используются весьма успешно. 0x00 - это Start, 0x08 - SyntaxCheck, 0x10 - NextChar, 0x18 - OutChar, 0x20 - CompareHLDE, и вот по 0x28 лежит FTestSign (узнать знак числа с плавающей точкой):

0028: 3A7201    FTestSign       LDA FACCUM+3    
002B: B7                        ORA A   
002C: C2DA09                    JNZ FTestSign_tail      
002F: C9                        RET

Снова вспомним: когда в регистрах по старшинству идут B (экспонента), C (знак и старшие биты мантиссы), D, E, то в памяти FACCUM это самый младший байт мантиссы, FACCUM+1 следующий, FACCUM+2 знак и старшие биты мантиссы, FACCUM+3 экспонента. И ещё напомним: это не стандартное представление числа с плавающей точкой "IEEE", его тогда ещё не было. Отличие - в расположении знака, он не в самом старшем бите, а сразу за экспонентой. И также тут нет денормализованных чисел, Inf и NaN.

Забежав немного вперёд, я выяснил: функция FTestSign возвращает в регистре A значение 0 (число строго ноль), 1 (число строго положительное) либо -1 (0xFF, если число строго отрицательное), и с соответствующими флагами (Z, если ноль в A, и S, sign, если A=-1).

При проверке знака первым делом проверяется байт экспоненты. Если там ноль (т.е и само число - НОЛЬ), то проверять сам знак не нужно, тем самым, "плюс ноль" и "минус ноль" (особенность всех ходовых представлений с плавающей точкой) оказываются "равны" между собой.

Если же в экспоненте не ноль (т.е не самое маленькое возможное значение), то вызывается "хвост" функции, FTestSign_tail:

09DA: 3A7101    FTestSign_tail  LDA FACCUM+2    
09DD: FE                        DB 0xFE 
09DE: 2F        InvSignToInt    CMA     
09DF: 17        SignToInt       RAL     
09E0: 9F                        SBB A   
09E1: C0                        RNZ     
09E2: 3C                        INR A   
09E3: C9                        RET

Первым делом загружаем байт FACCUM+2, где в старшем бите лежит знак, а потом мантисса (без старшей единицы, которая сидит "неявно").

Дальше идёт наложение инструкций. Дизассемблер предполагал, что управление передаётся на 0x09DE (функция InvSignToInt), но если мы прыгнули на FTestSign_tail, то при декодировании следующей инструкции 0xFE будет означать CPI (compare immediate), инструкция, занимающая 2 байта. Второй байт - число, с которым производится сравнение, в данном случае это будет 2F.

Сам результат сравнения нафиг не нужен, это лишь способ "перепрыгнуть" через ненужную инструкцию CMA (complement A - побитная инверсия регистра A).

Следующая команда: RAL, это циклический сдвиг влево через перенос (Carry). Так у нас знак попадает в Carry, а в регистр A вдвигается справа Carry, полученный как результат сравнения с 0x2F.

SBB A - ещё одна "странная инструкция": мы из A вычитаем A, но ещё вычитаем Carry, т.е "заимствование с прошлого разряда". Как результат: если знак был "+" (что кодируется нулём в старшем бите), в A установится ноль, иначе: 0xFF (он же 255, он же "-1").

RNZ - Return if Not Zero - ежели получилось "-1", то здесь же и возвращаемся (в нашем случае - в FMUL). Иначе ещё прибавляется единичка к A, т.е ответом становится "1".

Возвращаясь к FMUL, видим: если в "аккумуляторе" FACCUM лежал нуль, то просто завершаем процедуру (Return if Zero), ведь что ни умножай на ноль, так и останется ноль. В "современном представлении", это не совсем так: умножение нуля на NaN или Inf должно выдать NaN. Но здесь, в бейсике, подобной фигни не было.

Смотрим дальше:

08E7: 2E00              MVI L,00h       ;L=0 to signify exponent add
08E9: CD9B09            CALL FExponentAdd

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

099B: 78   FExponentAdd MOV A,B  
099C: B7                ORA A   
099D: CABA09            JZ FExponentAdd+31      
09A0: 7D                MOV A,L ;A=0 for add, FF for subtract.
09A1: 217201            LXI H,FACCUM+3  ;
09A4: AE                XRA M   ;XOR with FAccum's exponent.
09A5: 80                ADD B   ;Add exponents
09A6: 47                MOV B,A ;
09A7: 1F                RAR     ;Carry (after the add) into bit 7.
09A8: A8                XRA B   ;XOR with old bit 7.
09A9: 78                MOV A,B ;
09AA: F2B909            JP FExponentAdd+30      ;If
09AD: C680              ADI 0x80         
09AF: 77                MOV M,A  
09B0: CA2109            JZ PopHLandReturn        
09B3: CD370A            CALL FUnpackMantissas   
09B6: 77                MOV M,A 
09B7: 2B                DCX H   
09B8: C9                RET     
09B9: B7                ORA A   ;это строка FExponentAdd+30
; следующая строка - аккурат FExponentAdd+31
09BA: E1                POP H   ;Ignore return address so we'll end 
09BB: FAA408            JM Overflow     
09BE: AF        FZero   XRA A   
09BF: 327201            STA FACCUM+3    
09C2: C9                RET     

Загружаем в A экспоненту левой части и проверяем на ноль (т.е что в этот раз левая часть нулевая). Если так, перепрыгиваем на 0x09BA, где скидывается адрес возврата в процедуру FMUL (вслед за ним идёт адрес возврата туда, откуда был вызван FMUL). JM выполниться не может, поскольку когда результат ноль, флаг "минус" не установлен.

Следующие 3 строки уже использовались в FAdd / FSub, хотя мы забыли их включить в общее количество строк. Это установить нулевой результат. В этом стандарте чисел с плавающей точкой, это сделать легко: достаточно установить экспоненту в ноль. В более новых "IEEE", по-хорошему, надо вообще все байты установить нулевыми, чтобы не получилось денормализованных чисел.

Если же экспонента не нулевая, то всё-таки попробуем их сложить. Если бы экспонента была "без смещения", т.е "0" означало бы конкретно 2⁰, "1" - 2¹ и так далее, то можно было бы просто их сложить. Вот тута: https://en.wikipedia.org/wiki/Microsoft_Binary_Format
они, увы, сами себе противоречат. Вроде пишут о смещении в 128, дескать, значение 128 означает экспоненту "0", тогда 1 будет означать -127, а 255 будет означать +127. Но с другой стороны, иллюстрируют, что число "1" в этом формате записывается как 0x81 00 00 00, а вот "0,5" - уже 0x80 00 00 00. Видимо, здесь было принято записывать число в формате 2^e*0,1xxxxx, т.е "неявная единица" мантиссы - это не в точности единица, а вовсе даже 1/2. Что-то в этом есть: тогда после умножения двух мантисс результат по-прежнему будет лежать в диапазоне 0,25..1, что интерпретировать чуть проще. Диапазон представимых экспонент тогда - от -127 до +127. -128 в нашей записи означет "нулевое число".

Тем не менее, начинаем с простого сложения, ADD. Результат сохраняем в регистр B, а регистр A тут же крутим вправо через перенос. Тем самым, бит переноса оказывается в старшем бите A. Далее делается XOR с результатом сложения (сохранённым в B). Это очень похоже на сложение знаковых чисел, когда признаком переполнения становится совпадение бита переноса и старшего бита результата.

Попытаемся сообразить, как это работает. Настоящие экспоненты складываются так:

 (E1 - 128) + (E2 - 128) = (E1 + E2) - 256

При этом E1, E2 принимают значения от 1 до 255 (ноль означает нулевое число, и тогда мы сюда уже не попадаем). Следовательно, после сложения может получиться от -254 до +254. Но если получилось от -254 до -128, это не влезет в результат, поэтому надо будет записать ноль, а если от +128 до +254 - надо выругаться на переполнение.

в процессе умножения мантисс мы можем получить результат меньше ожидаемого: от 0,25 до 0,5, и тогда надо будет уменьшить экспоненту на единицу, а число подвинуть влево. Поэтому, теоретически, нам бы и +128 не стоило бы сразу сбрасывать со счетов!

Теперь переведём эти размышления в 8 бит результата и бит переноса. Оба будут нулевыми, если E1+E2 < 128, что соответствует экспоненте менее -128. Тут результатом XOR'а станет нулик в старшем разряде, что соответствует "положительному значению", поэтому JP FExponentAdd+30 (Jump if Positive) здесь исполнится. заметим: если E1+E2=128, то здесь мы ещё не прыгнем, хотя это соответствует непредставимой экспоненте

Оба бита будут единичными, если E1+E2 >= 384, что соответствует экспоненте большей или равной 128. Да, здесь мы тоже исполним прыжок. Иногда мы ошибочно заявим о переполнении, хотя на самом деле лишь ПРОМЕЖУТОЧНЫЙ результат получился из серии 2¹²⁸*0,01..., а после нормализации он бы вполне уместился. В общем, на самом краю динамического диапазона работать не стоит...

Прыгаем на один байт раньше, чем в случае нуля. Там строка ORA A устанавливает флаги в соответствии с нашим результатом сложения. Так что в первом рассмотренном случае результат будет "положительным", и мы закончим установкой нулевого результата. А во втором случае результат "отрицательный" - и мы прыгаем в Overflow, чтобы грязно выругаться!

Если же прыжка не было, то E1+E2 лежит в диапазоне от 128 до 383, и соответствует это экспонентам от -128 до +127, которые должны кодироваться значениями от 0 до 255. Т.е, чтобы получить правильный результат, нужно ещё отнять 128, или, что то же самое, прибавить 128. Команда ADI 0x80 делает именно это. Результат тотчас же помещается назад в "аккумулятор", по FACCUM+3. Но далее, если результат нулевой (что возможно лишь в одном случае, перед сложением у нас было 128 без переноса, что соответствует экспоненте -128), мы опять преждевременно заканчиваем выполнение, в этот раз через вызов PopHLandReturn. На регистры H/L нам плевать, главное, вышвырнуть из стека "лишний" адрес возврата, ведущий нас в FMUL, чтобы выполнение FMUL завершилось уже здесь.

FUnpackMantissas мы уже разбирали в FADD. Знак чисел заменяется на единицы, т.е "неявные единицы мантиссы" стали явными. Притом знак одного из чисел сохраняется в FTemp, а результат XNOR от двух знаков сейчас (по выходу из FUnpackMantissas) лежит в старшем бите A. Т.е, единица, если знаки совпадают, и нолик, если знаки различны. При этом регистры H/L в данный момент указывают на FTemp, где хранился знак первого операнда.

Здесь, только возвратившись из FUnpackMantissas, первым делом туда же в FTemp запихиваем A, поскольку нам при умножении чисел и нужно "знаки помножить". Сейчас там, выходит, будет лежать всегда противоположный знак.

DCX H - это чтобы (H,L) указывало не на FTemp, а на FAccum+3, т.е на экспоненту.

И - о чудо - НОРМАЛЬНЫЙ ЧЕЛОВЕЧЕСКИЙ RET! Всё, что идёт вслед за ним - мы уже рассмотрели, выход по "нулю" либо по переполнению.

Итак, мы снова в функции FMUL! Смотрим, что там дальше:

08EC: 79                MOV A,C 
08ED: 321709            STA FMulInnerLoop+13    
08F0: EB                XCHG    
08F1: 221209            SHLD FMulInnerLoop+8    
08F4: 010000            LXI B,0000h     
08F7: 50                MOV D,B 
08F8: 58                MOV E,B 
08F9: 215E08            LXI H,FNormalise+3      
08FC: E5                PUSH H  
08FD: 210509            LXI H,FMulOuterLoop     
0900: E5                PUSH H  
0901: E5                PUSH H  
0902: 216F01            LXI H,FACCUM

Тут у нас, кажись, пошёл самомодифицирующийся код!
В регистр A поместили старший байт мантиссы левого операнда (теперь уже с неявной единицей), а затем сохранили по адресу FMulInnerLoop+13, или 0x0917. Посмотрим на этот кусочек:

0916: CE00                  ACI 00  ;A=result mantissa high byte. This gets back to C

Регистр C уйдёт во второй байт, который пока нулевой. ACI - это Add Immediate with Carry, так что второй байт - тот самый Immediate, который прибавляется. Позже сообразим, для чего это надо.

Далее идёт XCHG, который меняет местами (D,E), где лежат наши младшие байты мантиссы и (H,L), где лежал указатель на FACCUM+3.

SHLD - это Store (H,L) Directly. То есть, регистр L отправляется по адресу FMulInnerLoop+8, или 0x0912, а регистр H - на единичку дальше, FMulInnerLoop+9, или 0x0913. Глянем, что там творится:

0911: 210000            LXI H,0000h

Ага, ещё одна Immediate инструкция. X здесь означает 16-битные операции (как впоследствии в x86 появились регистры AX, BX и тд, состоящие из AH и AL), в общем, это загрузить в (H,L) непосредственные значения, те самые, которые мы сюда запихиваем.

Далее в этом "подготовительном коде" мы зануляем регистры B,C,D,E. Благо, B уже не нужен (экспоненты мы уже сложили, результат лежит в FACCUM+3), а C,D,E мы уже запихали непосредственно внутрь программного кода. Полагаю, тут будет большущая площадка для складирования произведения двух мантисс.

В (H,L) помещаем адрес FNormalize+3, и тут же запихиваем его в стек. Это снова игры с tail call - мы добавили ещё один адрес возврата, чтобы при возвращении отсюда через Ret сначала попасть в FNormalize, а уже оттуда - вернуться туда, откуда вызывался FMUL.

Дальше ещё страньше: теперь мы ДВАЖДЫ помещаем на стек адрес FMulOuterLoop. И, наконец, в пару (H,L) загружаем FACCUM, это младший байт мантиссы "аккумулятора".

Похоже, мы подготовились к умножению... Вот оно:

0905: 7E  FMulOuterLoop MOV A,M ;A=FACCUM mantissa byte
0906: 23                INX H   ;
0907: E5                PUSH H  ;Preserve FACCUM ptr
0908: 2E08              MVI L,08h       ;8 bits to do
090A: 1F  FMulInnerLoop RAR     ;Test lowest bit of mantissa byte
090B: 67                MOV H,A ;Preserve mantissa byte
090C: 79                MOV A,C ;A=result mantissa's high byte
090D: D21909            JNC L0919       ;If that bit of multiplicand was 0, then skip over adding mantissas.
0910: E5                PUSH H  ;
0911: 210000            LXI H,0000h     ;
0914: 19                DAD D   ;
0915: D1                POP D   ;
0916: CE00              ACI 00  ;A=result mantissa high byte. This gets back to C
0918: EB                XCHG    ;in the call to FMantissaRtOnce+1.
0919: CDD708    L0919   CALL FMantissaRtOnce+1  
091C: 2D                DCR L   
091D: 7C                MOV A,H ;Restore mantissa byte and
091E: C20A09            JNZ FMulInnerLoop       ;jump back if L is not yet 0.
0921: E1 PopHLandReturn POP H   ;Restore FACCUM ptr
0922: C9                RET     ;Return to FMulOuterLoop, or if finished that then exit to FNormalise

Итак, загружаем сначала младший байт мантиссы правого операнда (FACCUM). Инкрементируем наш указатель, чтобы указывал на следующий по старшинству. Сохраняем его на стеке, поскольку сейчас применим L в качестве счётчика до 8.

Начинается внутренний цикл. Крутим регистр A вправо, через перенос. Младший бит остаётся в переносе. То, что осталось, сохраняем в H, а в A "заблаговременно" заносим значение из регистра C. Пока что он был занулён.

Если "текущий бит", очутившийся в переносе, был нулём, сложение на текущей итерации пропускаем.

Опять сохраняем регистры H,L, после чего загружаем в них два младших байта левого операнда, те самые, которые мы сохранили непосредственно здесь в коде.

Через DAD D мы к (D,E) прибавляем эти два байта, результат ложится в (H,L). Ранее сохранённые (H,L), где H это частично сдвинутый байт мантиссы, а L счётчик, сколько бит осталось, теперь ложатся в (D,E).

Ровно здесь мы к старшему байту результата, C, прибавляем старший байт левого операнда, сохранённый непосредственно в коде, ещё и с переносом из младших.

XCHG обменивает (D,E) и (H,L), так что "всё возвращается на круги своя", ну почти:

- в A лежит старший байт промежуточного результата,
- в D,E лежат следующие два байта промежуточного результата,
- в H лежит частично сдвинутый байт мантиссы второго числа,
- в L - сколько раз ещё нужно провести внутренний цикл.

На этом месте вызывается FMAntissaRtOnce+1. Эта штука у нас уже была, сдвиг мантиссы на единичку вправо. Но здесь мы, похоже, одну строчку пропускаем. Вот начало:

08D6: 79 FMantissaRtOnce        MOV A,C 
08D7: 1F                        RAR

Ну да, у нас уже старший байт в A сидит, а вот по окончании FMantissaRtOnce вернётся в C.

DCR L - вычитаем единичку из счётчика.

MOV A, H - готовимся к следующей итерации, возвращаем в A частично сдвинутый байт мантиссы.

JNZ FMulInnerLoop - заставляет нас провернуть 8 итераций по внутреннему циклу.

Да, картина начинает складываться. По сути, реализовано умножение "в столбик". Вроде такого:

(C,D,E) = xxxx_xxxx_FEDC_BA98_7654_3210
(FACCUM)= xxxx_xxxx_FEDC_BA98_7654_3210
----------------------------------------
          xxxx_xxxx_xxxx_xxxx_xxxx_xxxx

"Верхнее" число (оно же "левое") мы заблаговременно сохранили внутрь опкодов, поскольку оно будет всё время складываться "само с собой", только со сдвинутым. Соответственно, (C,D,E) обнулили.

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

А затем весь результат сдвигаем вправо на единичку, т.к у нас каждый раз выполняется лишь 24-битное сложение. Смещая результат вправо, мы "убираем из сумматора" те биты, к которым прибавлять уже нечего. По крайней мере, у нас есть ещё один вспомогательный регистр, куда вдвигаются самые младшие биты, что позволит потом довольно корректно провести округление.

Такой способ умножения, со сдвиганием результата вправо, наиболее экономичен, но не позволяет расширить умножение до "умножения с накоплением" (fused multiply-add). Ну здесь оно и не надо пока...

Когда заканчивается внутренний цикл, мы возвращаем из стека адрес очередного байта мантиссы, который мы возьмём из FAccum.

И здесь, наконец, стало ясно, почему FMulOuterLoop заложен в стек дважды. Так мы "возвращаемся" один раз, чтобы обработать второй байт, и ЕЩЁ один раз - обработать третий байт - а потом тем же Ret наконец-то попадаем в процедуру нормализации!

При этом последние две строчки кода, как видно, ещё и "по совместительству" являются процедурой PopHLandReturn, которую мы обсудили раньше.

На этом почти всё! Осталось посмотреть, куда мы а FNormalize прыгаем, мы же снова не в самое начало:

085B: DCB508 FNormalise CC FNegateInt   ; 
085E: 2600              MVI H,00h       ;
0860: 79                MOV A,C ;Test most-significant bit of mantissa
0861: B7                ORA A   ;and jump ahead if it's 1.
0862: FA7E08            JM FRoundUp     ;

Да, мы пропускаем первую команду, Call if Carry - она была нужна при вычитании мантисс, т.к знак мог поменяться, и тогда результату нужно было поменять знак. Но здесь такого произойти не может.

Здесь использовалась интерпретация мантиссы как числа от 0,5 до 1, но после умножения она может оказаться слишком маленькой, от 0,25 до 0,5, и тогда нужно сдвинуть её влево, а экспоненту уменьшить на единичку.

Ну а дальше ещё и округление, благо у нас результат умножения 32-битный, т.е сидят дополнительные 8 бит. Из них будет взят старший, и если он единичный - к 24-битной мантиссе ещё добавится единичка младшего разряда.


Этот код мне нравится. Видно, что и диапазон экспоненты здесь подбирали аккурат с тем расчётом, чтобы оно тут элегантно реализовывалось. Иногда эта штука пожалуется о переполнении там, где его быть не должно (переполняется лишь промежуточный, ненормализованый результат), так что можно заявить: реально диапазон представимых чисел не от 2^-128 до 2¹²⁷, а "всего лишь" до 2¹²⁶. Но думаю, это пережить можно.

Оставить комментарий

Популярные посты:

• Ранее
• Архив