Изучаем FDIV для Altair Basic

nabbla1 — 04.03.2025 Сложение/вычитание смотрели, умножение смотрели, теперь вот до деления дошли. Традиционно деление считается самой сложной арифметической операцией, вспомнить хотя бы ошибку FDIV в пентиумах.

В коде для Altair Basic (под процессор Intel 8080) сама процедура деления с плавающей точкой занимает 107 байт, это не считая используемых не только делением процедур, вроде определения знака числа, сложения/вычитания экспонент, распаковки мантисс. Как и для умножения, вовсю используется самомодифицирующийся код, есть и наложение команд в одном месте, и выход "сквозь" другую процедуру, т.е CALL FRoundUp, RET заменяется на JMP FRoundUp.

Код, очевидно, рабочий, но в сравнении с умножением (FMUL) граничные случаи рассмотрены ещё более "спустя рукава". Слишком большие числа могут "преждевременно" вызвать переполнение (хотя ответ всё ещё был представим), слишком маленькие - обратиться в ноль (хотя какая-то точность ещё могла сохраниться), но что самое нехорошее - этот код иногда путает большие и маленькие числа! Там, где должно было выйти переполнение, вдруг получится нулевой результат, а где должна была произойти полная потеря точности (до нуля) - вдруг выскакивает переполнение!

Попробуем разобраться...


Листинг берём отсюда: https://github.com/option8/Altair-BASIC/blob/master/BASIC%20disassembly-source.lst
деление начинается по адресу 0x92F.

Но прежде, чем рассматривать сам код, расскажем словами, что там происходит, на примере 2.0 / 3.0.

Число 2.0 представляется как экспонента (первый байт) 130 = 0x82 и мантисса (следующие 3 байта): 0x00 00 00. В старшем бите лежит знак, равный 0 (то есть "плюс"), но когда дело дойдёт до деления мантисс, на это место встанет "неявная единица". Т.е получается число 2^130-128 × 0.10000...₂.

Число 3.0 представляется экспонентой 130 = 0x82 и мантиссой 0x40 00 00. Когда дело доходит до деления, на место старшего бита встаёт "неявная единица" (становится 0xC0 00 00), т.е получается число 2^130-128 × 0.110000....₂.

Первым делом производится вычитание экспонент, но для этого используется уже применявшийся в FMUL код FExponentAdd. Но в регистр L заблаговременно заносится значение 0xFF (оно же 255, оно же все единички), чтобы провести XOR этого числа с экспонентой делителя. Получается инверсия: 0x82 = 1000 0010₂ превращается в 0111 1101₂ = 0x7D.

Далее прибавляется 0x82 (экспонента делимого) и 0x7D (инвертированная экспонента делителя), что даёт 0xFF = 255, без флага переноса. Различие бита переноса и старшего бита результата говорит: переполнения не произошло, всё нормально. Далее к результату прибавляется 0x80 = 128, что даёт 127 + бит переноса.

Но затем, уже по выходу из FExponentAdd к экспоненте ещё прибавляется 2, что даёт 129 = 0x81. Т.е пока что предполагается, что результат деления будет лежать в диапазоне 1..2.

После этого приступаем к делению мантисс. Мантисса делителя, 0xC0 00 00, целиком помещается внутрь исполняемого кода, в immediate инструкции SUI и SBI (Sub Immediate и Sub Immediate with Borrow). Делимое расширяется на один байт влево, изначально он нулевой:
0x00 80 00 00. Наконец, на частное выделяется 3 байта, изначально они инициализируются нулями.

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

 00 80 00 00
-00 C0 00 00

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

Делимое сдвигается влево, превращаясь в 0x01 00 00 00.

В конце итерации проверяется, равно ли частное в точности нулю. Так и есть: значит, из экспоненты вычитаем единичку, теперь это 128=0x80 (т.е результат лежит в диапазоне 0,5..1, так и есть).

Начинаем следующую итерацию. Снова из текущего делимого вычитаем делитель:

 01 00 00 00
-00 C0 00 00
-------------
 00 40 00 00

Результат положительный, следовательно, в частное вдвигается единичка, получается 00 00 01. Делимое сдвигается влево, давая 00 80 00 00. Частное уже не ноль (и заведомо нулём не будет), так что экспоненту больше не трогаем.

Дальше, как можно заметить, у нас всё "зацикливается" - делимое всё время скачет между 00 80 00 00 (не получается вычесть делитель) и 01 00 00 00 (вычитается, даёт 00 40 00 00, которое снова двигается влево и снова получается 00 80 00 00), при этом в частное вдвигается последовательность 0, 1, 0, 1 и так далее.

Так продолжается, пока у нас не установится в единицу самый старший бит частного. Т.е оно принимает вид 1010 1010 1010 1010 1010 1010, или 0xAA AA AA. Как только это произошло - мы возвращаемся из FDIV, но не сразу, а через процедуру FRoundUp, причём в качестве младшего бита, определяющего, нужно ли к ответу прибавить единичку (округлить вверх), выступает бит переноса с последней итерации. Там единичка, поэтому финальная мантисса: 0xAA AA AB, или 11 184 811. Поделив её на 256³ = 0x01 00 00 00, получаем 0,666666686534881591796875.
Без округления "вверх" получилось бы 0,66666662693023681640625, что чуть дальше от истины.

Вот как-то так оно и работает. Теперь чуть подробнее о коде, и что мне там не нравится.

Вот самое начало, по адресу 0x92F:

092F: C1        FDiv    POP B    
0930: D1                POP D   
0931: EF                RST 05  ; FTestSign      
0932: CAD301            JZ DivideByZero 
0935: 2EFF              MVI L,0xFF       
0937: CD9B09            CALL FExponentAdd       
093A: 34                INR M    
093B: 34                INR M

"вступительный код" - здесь мы вспоминаем, что правый операнд (в данном случае делитель) лежит в "аккумуляторе" FACCUM, располагаемом в памяти, а левый операнд (делимое) - пока что на вершине стека. Первым же делом он оттуда извлекается в регистры B,C,D,E (в B экспонента, в C знак и старшие биты мантиссы, E - самые младшие биты мантиссы).

Тут же с помощью "прерывания" (RST, Restart) вызывается проверка знака числа в аккумуляторе и по совместительству проверка на ноль. Ежели число оказалось нулевым - прыгаем на DivideByZero, где ругаемся на пользователя.

В противном случае в регистр L заносим 0xFF, оно же 255, оно же "все единицы", и вызываем FExponentAdd. Подробно он разбирался в прошлой части, про FMUL. Осуществляется XOR экспоненты делителя с регистром L, что равносильно смене знака и вычитанию единички. К примеру, 0 превратится в 255 и будет по своему действию эквивалентно "-1".

Если экспонента делимого - ноль (т.е и сам делитель - ноль, так уж условились) - прямо там, внутри FExponentAdd будет установлен нулевой результат и произведён возврат туда, откуда была вызвана FDiv, т.е на этом деление будет завершено.

В противном случае из экспоненты делимого вычитается экспонента делителя (и ещё минус 1). Самое время снова громко подумать... Мы пришли к выводу, что в Altair Basic числа представляются следующим образом: 2^E-128×0.1mmmm.... (E-записанная в первом байте экспонента, 0..255, 0 особый случай, mmmm - биты мантиссы, 23 штуки).
Т.е экспонента записана со смещением 128, при этом мантисса со вписанной в неё "неявной единицей" может принимать значения от 0,5 почти до 1. Так оно оказывалось удобнее при умножении: обычное "умножение в столбик" можно было интерпретировать в том же диапазоне, "до 1", разве что если там слишком маленькое значение, от 0,25 до 0,5 - производился сдвиг влево с вычитанием единички из экспоненты.

В этот раз, когда делим две мантиссы, принимающие значение от 0,5 до 1, то ответ может лежать от 0,5 до 2. При этом экспоненты вычтутся вот так:

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

Что и логично: при вычитании наше "смещение" ушло. Когда E1, E2 принимают значения от 1 до 255 (ноль означает что число нулевое - эти случае мы уже отсекли), итоговая экспонента может равняться от -254 до +254. Тут всё в точности как при умножении, как это ни странно: значения от -254 до -128 у нас окажутся непредставимы - надо будет объявить число нулём. А значения от +128 до +254 объявить переполнением. (это без учёта мантиссы - по-хорошему и её бы проверить, но здесь всё-таки позволили себе немного снизить динамический диапазон ради упрощения и убыстрения кода)

В действительности происходит немного не так: считается скорее E1-E2-1, поскольку E2 лишь "инвертировали", а истинная смена знака требует ещё единичку прибавить.

Рассмотрим 4 случая:
1. переноса не было, старший бит результата нолик, т.е получилось число от 0 до 127. Такое может случиться, если E1 изначально был от 1 до 127, а E2 был от 128 до 255. Например, E1=127 (соотв. экспоненте "-1"), E2=255 (соотв. экспоненте "+127"). Если честно вычесть, получится -128 - это уже непредставимо. Но у нас E2 проинвертируется (получится 0), и затем к 127 прибавится ноль, так и останется 127. Отсутствие переноса в данном случае означает, что результат на 256 меньше, т.е -129.

Заметим, "настоящая" экспонента в этом случае - от -128 до -254, то есть на этот раз вроде как всё хорошо отловили.

2. переноса не было, старший бит результата единичка, т.е получилось число от 128 до 255. Снова, отсутствие переноса в данной ситуации - это неявное вычитание 256 из результата, т.е на деле у нас как будто получились значения -128..-1 (а "реальные экспоненты": -127..0), которые вполне представимы. Например, E1=E2=128 (экспонента "0"). Инвертируем E2, оно превращается в 127. Складываем их - получается 255 без переноса, что надо интерпретировать как "-1". Нормальный исход, можно продолжить деление.

3. Был перенос, старший бит нолик, т.е число от 0 до 127. Это и есть экспоненты от 0 до 127, но если ещё точнее - от 1 до 128, поскольку мы там ненароком единичку вычли. Тоже пока сойдёт, можно продолжить деление.

4. Был перенос и старший бит единица, т.е число от 128 до 255. Это ПЕРЕПОЛНЕНИЕ, что-то из серии E1=255, а E2=1. Проинвертировали E2 - получилось 254. Сложили их - получили 253 и флаг переноса. И действительно, в этой ситуации будет выдано переполнение.

В общем, как это ни странно, та же самая логика, которую применяли при умножении чисел, подходит здесь при делении.

Потом ещё прибавляется 128, и если вдруг получился ноль, назначим результат нулём и прекратим деление. Да: так бывает, если переноса не было и получилось число 128, которое соответствовало настоящей экспоненте в -127, вполне себе представимой. Но увы, "для единообразия" получается потеря точности - число ещё было представимо, но стало нулём.

Вернувшись из FAddExponents, мы ещё дважды прибавляем к получившейся экспоненте единичку. Одну "по делу", вторую, видимо, чтобы интерпретировать получающуюся мантиссу как принимающую значения 0,25..1, а не 0,5..2, как происходит, когда число 0,5..1 делишь на другое число 0,5..1.

К этому прибавлению двойки у меня наибольшие претензии! Т.к тут вдруг нет проверок на переполнение. То есть, внутри FExponentAdd у нас мог получиться результат 254 или 255, а здесь, после прибавления двойки - оно превратится в 0 или 1, и МОЛЧОК. Теоретически, это куда страшнее, чем ошибка в Pentium FDIV, но там где "железячники" попали на полмилииарда убитых енотов, с мелкого и мягкого как с гуся вода. Ну да, ничего жизненно важного на бейсике писаться не должно было, так что может и есть в этом сермяжная правда.

Поехали дальше. Подготовительные операции перед циклом:

093C: 2B                DCX H   
093D: 7E                MOV A,M 
093E: 326009            STA L095F+1     
0941: 2B                DCX H   
0942: 7E                MOV A,M 
0943: 325C09            STA L095F-3     
0946: 2B                DCX H   
0947: 7E                MOV A,M 
0948: 325809            STA L095F-7     
094B: 41                MOV B,C 
094C: EB                XCHG    
094D: AF                XRA A   
094E: 4F                MOV C,A 
094F: 57                MOV D,A 
0950: 5F                MOV E,A 
0951: 326309            STA L095F+4

Пара (H,L) ранее указывала на экспоненту, теперь вот вычли единичку - и теперь указывает на старший байт мантиссы.
Загрузили оттуда значение и разместили "внутри исполняемого кода", в команду SBI 00h. Т.е вместо 00h здесь будет лежать конкретный байт делителя.

Потом ещё раз вычитаем 1 из указателя (H,L), загружаем средний байт мантиссы делителя - и размещаем его в другую команду SBI 00h.

Третий раз то же самое - теперь это младший байт, его размещаем в команду SUI 00h.

Делимое лежало в регистрах C (старший байт мантиссы), D (средний), E (младший). В регистре B лежала экспонента, но она уже не нужна (её мы посчитали почти что, результат лежит в FACCUM+3). Теперь старший байт отправляется в регистр B, а средний и младший (D,E) с помощью XCHG помещаются в (H,L).

XRA A - попросту ОБНУЛЕНИЕ регистра А. И далее через него обнуляются C,D,E - теперь там будет лежать частное, изначально нолик.

И наконец, ещё один "самомодифицирующийся" код - обнуляем "старший байт делимого". Он также лежит прямо внутри исполняемого кода - регистров катастрофически не хватало!

Далее начинается сам цикл деления:

0954: E5      FDivLoop  PUSH H  
0955: C5                PUSH B  
0956: 7D                MOV A,L 
0957: D600              SUI 00h ;здесь вместо 00 будет младший байт делителя
0959: 6F                MOV L,A 
095A: 7C                MOV A,H 
095B: DE00              SBI 00  ;здесь вместо 00 будет средний байт делителя
095D: 67                MOV H,A 
095E: 78                MOV A,B 
095F: DE00      L095F   SBI 00  ;здесь вместо 00 будет старший байт делителя
0961: 47                MOV B,A 
0962: 3E00              MVI A,00h ;здесь вместо 00 будет лежать САМЫЙ старший байт делимого, непрерывно обновляемый     
0964: DE00              SBI 00  
0966: 3F                CMC     
0967: D27109            JNC L0971       
096A: 326309            STA L095F+4h    
096D: F1                POP PSW 
096E: F1                POP PSW 
096F: 37                STC     
0970: D2                DB 0xD2 ;JNC ....       
0971: C1        L0971   POP B   
0972: E1                POP H   
0973: 79                MOV A,C 
0974: 3C                INR A   
0975: 3D                DCR A   
0976: 1F                RAR     
0977: FA7F08            JM FRoundUp+1   
097A: 17                RAL     
097B: CD9008            CALL FMantissaLeft      
097E: 29                DAD H   
097F: 78                MOV A,B 
0980: 17                RAL     
0981: 47                MOV B,A 
0982: 3A6309            LDA L095F+4h    
0985: 17                RAL     
0986: 326309            STA L095F+4h    
0989: 79                MOV A,C 
098A: B2                ORA D   
098B: B3                ORA E   
098C: C25409            JNZ FDivLoop    
098F: E5                PUSH H  
0990: 217201            LXI H,FACCUM+3  
0993: 35                DCR M   
0994: E1                POP H   
0995: C25409            JNZ FDivLoop    
0998: C3A408            JMP Overflow

Делимое (регистры C, H, L) "на всякий случай" сохраняется в стек.

Далее начинается классическое вычитание "в столбик", справа налево. Поскольку делитель занимает 3 байта, а делимое мы расширили до 4 байт, то в старшем байте используется просто SBI 0, что означает "вычесть единичку, если было заимствование".

CMC - это CompleMent Carry, т.е инвертировать бит переноса. Дело в том, что по результатам этого SBI 0 будет C=0, если "результат положительный", и C=1, если "результат отрицательный". А в частное хочется запихнуть строго наоборот - единичку, если вычитание было "успешно" и нолик, если нет.

Если результат получился меньше нуля, мы прыгаем на метку L0971, где восстанавливаем делимое, точнее его 3 байта, C,D,E. Самый старший байт восстанавливать и не надо, т.к его мы назад не записывали.

Если же прыжка не было, то как раз записываем самый старший байт делимого, а затем "выкидываем" сохранённые результаты из стека, но попросту "выкинуть" их тут нельзя (а через декремент SP тянуть очень уж муторно, на 4 позиции), поэтому они выпихиваются в сколько-нибудь "ненужные" регистры A и PSW (Program State Word, т.е как раз все флажки). Единственное, что при этом мы "теряем" флаг переноса, поэтому снова ручками устанавливаем его в единицу с помощью STC (SeT Carry).

И затем идёт наложение команд: чтобы "перпрыгнуть" через две команды POP, они превращаются в адрес для прыжка JNC, но только мы только что установили Carry=1, так что прыжка не будет, и мы успешно приступим к исполнению кода после двух POP.

Тут мы проверяем, достаточно ли итераций уже прошло. Для этого нужно определить старший бит частного. Если он уже единичка, значит - пора завершать. Старший байт частного - регистр C. Загружаем его в A, а затем прибавляем и вычитаем единичку, тупо, чтобы установились флажки.

Команда RAR загоняет флаг переноса в старший бит A. Это нужно внутри процедуры FRoundUp: с её точки зрения, в A лежал "дополнительый" (самый младший) байт мантиссы при сложении/вычитании, используемый для сколько-нибудь корректного округления.

Следующая команда, JM FRoundUp - это и есть основной выход из процедуры FDIV! Под "минусом" имеется в виду единица в старшем бите частного, означающим: провели сколько надо итераций. А уже внутри FRoundUp произойдёт округление, а затем возвращение знака числа на своё место, и оттуда уже вернёмся туда, откуда вызывался FDIV.

Если же мы ещё не закончили - "отменяем" действие RAR при помощи RAL (циклический сдвиг влево через перенос).
Далее частное сдвигается влево при помощи процедуры FMantissaLeft, причём справа "вдвигается" только что полученный бит результата, лежащий у нас во флаге переноса.

Далее сдвигается влево всё делимое. Два младших байта сдвигаются с помощью DAD H, т.е к (H,L) прибавить (H,L). Следующий байт лежит в регистре B - его пихаем в A и производим циклический сдвиг влево через перенос. Результат возвращаем в B, и затем то же самое проделываем с самым старшим байтом, лежащим внутри исполняемого кода, в команде MVI xx, т.е MOVe Immediate.

Проверяется равенство частного нулю, путём объединения всех трёх байтов через OR. Только если они все нули - результатом станет ноль. Тогда прыжок JNZ FDivLoop не сработает, вместо этого сохраним 2 байта делимого (H,L) в стек, а сами пока используем (H,L) как указатель на экспоненту результата, FACCUM+3. Вычтем единичку из экспоненты. Вернём на место делимое из стека, и затем, если экспонента не обратилась в ноль, начинаем очередную итерацию цикла.

Меня удивляет последняя строчка, JMP Overflow Какой же это, нахрен, Overflow, если, напротив, тут число стало СЛИШКОМ МАЛЕНЬКИМ???

Притом, поскольку экспонента и так стала нулём, это число уже по определению НОЛЬ. Т.е вместо JMP Overflow можно было просто написать Ret - и это было бы более логичным поведением, ещё и меньше места заняло бы, т.к RET занимает всего один байт.

Не сказать, что на эти особенности работы (обнуление вместо переполнения или наоборот, переполнение вместо обнуления) так уж легко наткнёшься, когда пишешь программки на бейсике. Но всё же наткнуться на них куда проще, чем на ошибку Pentium FDIV. Достаточно взять непрерывно растущее или непрерывно убывающее число - и обе ошибки проявятся ясно как день. Например, захочется школьнику посчитать какой-нибудь ряд, может число Пи попытаться найти или экспоненту какую. И очередное слагаемое ВНЕЗАПНО выдаст переполнение, хотя казалось бы, они стремительно убывают, в знаменателе факториал, как-никак!


Осталось рассмотреть квадратный корень и синус, в этом листинге это все математические функции, что есть.

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

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

• Архив