Изучаем квадратный корень для Altair Basic

топ 100 блогов nabbla105.03.2025

❝ Да, я взбунтовался! Я взбунтовался! И вы вовсе не величайший из королей, а всего лишь выдающийся, да и только! Что, съел? Выдающийся, да и только! И вы вовсе не по заслугам именуетесь почётным святым. Вы отшельник, подвижник, но не святой. Не святой! Нет! ❞


C четырьмя арифметическими операциями (сложение/вычитание, умножение, деление) разобрались, теперь ещё квадратный корень рассмотрим. В Basic он традиционно называется SQR (SQuare Root).

Эта функция относительно короткая, 61 байт, используется метод Ньютона. Судя по всему, считает точно (насколько позволяет одинарная точность), но по непонятной причине умножение на 0,5 здесь реализовано "в лоб", для чего "медленно и печально" набирается 4-байтовое число и вызывается FMUL. Декремент экспоненты вышел бы и существенно быстрее, и даже компактнее.

Рассмотрим работу на примере "корня из двух", а потом чуть подробнее изучим код. В этот раз совсем "низкоуровневых" фишек не наблюдается.


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

Допустим, мы хотим найти корень из двух. Двойка в формате чисел Altair Basic представляется как 0x82 00 00 00. Т.е экспонента 130, но там используется смещение 128, т.е в действительности только 2. Мантисса имеет неявную единицу сразу после запятой, а за ней в нашем случае идёт 23 нуля. Получается 22 × 0.12 = 2.

Первым делом проверяется знак числа - он положительный, и число не нуль. Значит, надо продолжать.

Затем экспонента делится на два с помощью циклического сдвига вправо через перенос (в переносе заведомо ничего нет - так получается, когда проверка знака завершается с положительным знаком). Получаем 0x41, а в бит переноса вдвигается нолик - младший бит экспоненты.

Это значение, 0x41, сохраняем на будущее, а пока что берём новое значение, 0x40, и сдвигаем его влево через перенос. Тем самым, задвинули сюда младший бит экспоненты (в нашем случае нолик) и получилось 0x80, что означает "от 0,5 до 1". Если бы младший бит был единицей, вышло бы 0x81, или "от 1 до 2".

Эта новая экспонента, 0x80, вставляется на место старой в исходное число, превращая двойку в 0,5. Т.е чётные степени двойки в экспоненте "вынесены за скобки", а к этому моменту число будет точно лежать в диапазоне от 0,5 до 2, так что большого количества итераций не требуется.

В качестве начальной итерации для метода Ньютона берётся само же число, в нашем случае 0,5.

Напомним, при начальном приближении x0, когда берём квадратный корень из числа A, надо выполнять итерации:

xn+1 = (xn + A/xn) / 2

Итераций фиксированное количество: ЧЕТЫРЕ.

Первая итерация довольно дурацкая: результатом становится (A+1)/2, в нашем случае 0,5 ("нулевое приближение") превращается в 0,75.

На второй итерации получается 0,708333...
На третьей итерации 0,7071078, и, наконец, на четвёртой - 0,7071068. Это правильное значение для корня из 1/2, насколько хватает разрядности.

Остаётся лишь "домножить" этот ответ на вынесенные за скобки степени двойки. Для этого возвращаем из стека значение 0x41, прибавляем к нему текущую экспоненту 0x80 (получается 0xC1). И наконец, прибавляем 0xC0 (фиксированное значение), что даёт 0x81 и никому не нужный перенос. Т.е значение "домножилось на два" (перешли от 0x80 к 0x81), так что ответом становится 1,414214. Это правильный ответ, насколько хватает разрядности.

Теперь всё-таки посмотрим на код. Значение, от которого надо взять корень, лежит в FACCUM. Сначала проверки на знак и на нулевое значение:
0C21: EF        Sqr     RST 05  ; FTestSign     ;
0C22: FA9804            JM FunctionCallError;
0C25: C8                RZ      ;


Затем танцы с бубном вокруг экспоненты:
0C26: 217201            LXI H,FACCUM+3  ;
0C29: 7E                MOV A,M ;
0C2A: 1F                RAR     ;
0C2B: F5                PUSH PSW        ;
0C2C: E5                PUSH H  ;
0C2D: 3E40              MVI A,40h       ;
0C2F: 17                RAL     ;
0C30: 77                MOV M,A ;


Напомним: PUSH PSW записывает и флаги, и регистр A, в данном случае нужен именно он. Регистры (H,L) сохраняем из экономии одного байта: без это надо было бы напрямую загрузить адрес (как первая команда в этом блоке, LXI H, FACCUM+3, занимает 3 байта), а так один байт сохранить и ещё один байт вовремя извлечь.

По окончании этих строк, в стеке лежит почти готовое финальное значение экспоненты, а исходное число (в FACCUM) теперь принимает значения от 0,5 до 2.

Ещё подготовка к циклу:
0C31: 217401            LXI H,FBUFFER   ;
0C34: CD290A            CALL FCopyToMem ;
0C37: 3E04              MVI A,04h       ;


В FBUFFER будет храниться неизменное значение, "из которого надо извлечь квадратный корень". А оно же, только в FACCUM - это "нулевое приближение", в дальнейшем оно изменится.

И наконец, в этом месте в A счётчик итераций, 4 штуки.

Рассмотрим цикл:
0C39: F5        SqrLoop PUSH PSW        ;
0C3A: CD020A            CALL FPush      ;
0C3D: 217401            LXI H,FBUFFER   ;
0C40: CD200A            CALL FLoadBCDEfromMem   
0C43: CD3109            CALL FDiv+2     
0C46: C1                POP B   
0C47: D1                POP D   
0C48: CD1208            CALL FAdd+2     
0C4B: 010080            LXI B,8000h     
0C4E: 51                MOV D,C 
0C4F: 59                MOV E,C 
0C50: CDE508            CALL FMul+2     
0C53: F1                POP PSW 
0C54: 3D                DCR A   
0C55: C2390C            JNZ SqrLoop     


Особенность почти всех здешних процедур для работы с плавающей точкой - они затирают все регистры. Так что понимаем: счётчик в регистре A не имеет ни единого шанса уцелеть - упихиваем его в стек.

Затем сохраняем в стек текущее приближение к квадратному корню.

Снова устанавливаем (H,L), чтобы указывали на FBUFFER (там хранится число, из которого извлекаем корень), и загружаем его в регистры B,C,D,E.
И затем делим его на текущее приближение. Тут повсюду используются FDIV+2, FADD+2, FMUL+2, поскольку первые 2 байта - вытаскивание из стека. А нам этого не надо, у нас уже всё в BCDE лежит.

В FACCUM теперь лежит A/xn, а в BCDE загружаем xn, и затем складываем их.

Теперь уже в FACCUM лежит числитель нашего выражения, остаётся поделить его на 2. Здесь это выполнено как "честное" умножение на число с плавающей точкой "0,5". В регистры (B,C) загружается значение 0x8000, а затем нолик из регистра C также копируется в D, E. Так получается число 0x80 00 00 00, т.е та самая "0,5". И после этого вызывается FMul+2, т.е умножение, когда второй множитель не надо извлекать из стека, он уже сидит в B,C,D,E.

Всего, чтобы записать это "умножение на 0,5", нам понадобилось 8 байт. Можно было бы сделать проще:

;осторожно, этот код я не проверял на эмуляторе, может что-то и напутал!
DCX H
DCR M


Поясняю: вот заключительное место при сложении с плавающей точкой, "основной выход":
0886: 46                MOV B,M ;B=exponent
0887: 23                INX H   ;
0888: 7E                MOV A,M ;A=FTEMP_SIGN
0889: E680              ANI 0x80        ;
088B: A9                XRA C   ;Bit 7 of C is always 1. Thi
088C: 4F                MOV C,A ;
088D: C3120A            JMP FLoadFromBCDE       ;Exit via copying BCDE to FACCUM.


Последнее, что происходит - в число закладывается правильный знак из FTEMP_SIGN, этот байт лежит сразу после FACCUM+3, где лежит экспонента. Далее выход производится через FLoadFromBCDE, который оставляет пару (H,L) неприкосновенной.

Так что нам остаётся уменьшить (H,L) на единичку, чтобы эта пара указывала на FACCUM+3 (на экспоненту нашего приближения), и затем уже из экспоненты вычесть единичку. ДВА БАЙТА ВМЕСТО ВОСЬМИ, и ускорение работы на десятки процентов, т.к умножение, хочешь не хочешь, 24 итерации, вот ни разу не коротких, займёт!

Тут можно не бояться всевозможных граничных ситуаций, т.к диапазон чисел, с которыми мы тут работаем, очень ограничен, от 0,5 до 2. Оно и в ноль ни в жисть не обратится, и переполнения не вызовет, и экспонента не может ВНЕЗАПНО уменьшиться до нуля, т.к принимает очень узкий диапазон значений: 0x80 либо 0x81.

Почему так не сделали - понятия не имею. Можно было бы на "освободившемся месте" сделать более качественное "нулевое приближение", что-то на манер FastInvSqrt, и сэкономить одну итерацию, хотя здесь из-за того, что в здешнем формате чисел с плавающей точкой между экспонентой и мантиссой затесался знак, это выходит далеко не так красиво...

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

Наконец, заключительные операции:

0C58: E1                POP H   
0C59: F1                POP PSW 
0C5A: C6C0              ADI 0xC0        
0C5C: 86                ADD M   
0C5D: 77                MOV M,A 
0C5E: C9                RET     


Возвращаем паре (H,L) адрес FACCUM+3 (там лежит экспонента нашего числа). А POP PSW загружает в регистр A старое значение, сдвинутое вправо. Поясню: у нас экспонента записывается со сдвигом 128. Поэтому если сидело значение E, то реальная экспонента E-128. Когда извлекается квадратный корень. эта "реальная экспонента" должна поделиться пополам, но записаться опять со смещением, а именно:

(E-128)/2 + 128 = E/2 + 64.

Пополам мы поделили. Тут бы прибавить эти 64 = 0x40, да и дело с концом. Но ещё младший бит потерять нельзя, поэтому прибавляем "текущую" экспоненту (128 или 129, т.е 0x80 или 0x81), но теперь ещё этот 0x80 "лишний" взялся, надо его "нейтрализовать", прибавив ещё 0x80 к нашему 0x40. Вот и получается 0xC0.

И в кои-то веки нормальный человеческий RET, иногда в этом коде прямо забываешь, что такое случается, без Tail call и условных возвратов (одна из фишек 8080, которая позже исчезла).

Если кто считает, что я придираюсь, а может самоутверждаюсь таким образом - см. эпиграф. Всё-таки этот код был написан совсем небольшим коллективом и за очень ограниченное время (хоть и не за 2 недели, у каждого участника уже были свои наработки, которые здесь пошли в ход. Это как собака Лайка, полетевшая в космос 3 ноября 1957 года после первого спутника 4 октября 1957, потому как подобные капсулы для собак уже пускали в геофизических ракетах много лет кряду, что не умаляет бешеного темпа, с которым тогда развивалась космонавтика). И хотя сейчас очень легко впечатлиться ВЕЛИКИМИ ДРЕВНИМИ, вместившими целый интерпретатор Бейсика в 4 килобайта, вместе со всей математикой (на фоне нынешних сред разработки, идущих в ГИГАБАЙТЫ), полезно всё ж помнить - они тоже были людьми из плоти и крови. Эти знания никуда не ушли, наоборот, приумножились, так что можно сделать и чуточку лучше, было бы желание!

Вообще, мне кажется, по ряду причин нас в скором времени ожидает очень основательная ревизия всего "стека компьютерных технологий", ибо вниз уже некуда (тёмный кремний пришёл, флэшки становятся всё более и более "аналоговыми", то бишь с множеством уровней, считываемых через АЦП, жёсткие диски записывая одну дорожку портят соседние и обязаны их "корректировать" - явно уже выскребаем остатки), так что экстенсивное развитие сменится интенсивным - зная, что нам по итогу нужно от этого хлама, сделать "начисто", как положено!


На днях расскажу о "быстром извлечении квадратного корня". На манер "обратного квадратного корня" из Quake, но даже проще! Один сдвиг вправо на единичку и одно сложение 32-битных чисел позволяет получить ответ с точностью не хуже 3%. А по математике Altair Basic осталось ещё синус рассмотреть. В 4 КБ версии, которую изучаю, похоже, других функций нет.

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

Предыдущие записи блогера :
Архив записей в блогах:
Октябрь 2022г. Нижний Новгород. Нижегородский кремль. Северная (Ильинская) башня. Координаты: 56°19'43"N 43°59'53"E ...
СКИДКА 25% на весь раздел для волос Хотела взять голубые спиральки к джинсам, да прощёлкала клювом. Но нежно-принцессочные (бело-перламутровые и розово-нюдовые) ещё пока держатся и доступны к покупке. Они нежно-девичьи, прямо-таки принцессочные. Очень и очень красивые. ...
Одно к одному? Но мне больше всего понравилось, какая смешная глушилка появилась в промо сразу после публикации новости о разводе. Теперь тому, кто рвется в ТОП с горячей новостью, придется раскошелиться:) ...
Экономисты Института исследований и экспертизы ВЭБ считают, что для роста численности населения недостаточно увеличить расходы на здравоохранение. По их мнению, власти РФ должны также оказать «меры поддержки трудовой и образовательной миграции», пишут «Ведомости». Изменить ситуацию, ...
в последнее время меня подключили уйма народа) давайте познакомимся, что-ли)))))) Эт я Алена Хозина и мой кот Рыжий) так, что фотки-в ...

A PHP Error was encountered

Severity: Warning

Message: mysqli::mysqli(): (08004/1040): Too many connections

Filename: libraries/Output.php

Line Number: 246

A PHP Error was encountered

Severity: Warning

Message: mysqli::query(): Couldn't fetch mysqli

Filename: libraries/Output.php

Line Number: 250