Квадратный корень на i8080

[shoorick]

Итерации по Ньютону. Может, кому пригодится ;)


;-----------------------------------------------------------------------
; HL = SQRT(HL)
;-----------------------------------------------------------------------
sqrt:
;-----------------------------------------------------------------------
mov a,h
ana l
inr a
jz .FFFF
;-----------------------------------------------------------------------
xchg
.loop:
mov b,h
mov c,l
push d
push b
lxi h,0
call gg_d32a
pop h
push h
dad d
mov a,h
rar
mov h,a
mov a,l
rar
mov l,a
;-----------------------------------------------------------------------
pop b
pop d
mov a,b
cmp h
jnz .loop
mov a,l
cmp c
jnz .loop
;-----------------------------------------------------------------------
ret
;-----------------------------------------------------------------------
.FFFF:
mov h,a
ret
;-----------------------------------------------------------------------


Использует подпрограмму деления Д32А из книжки Гуртовцева и Гудыменко


;-----------------------------------------------------------------------
; HLDE / BC = DE (HL)
; CY=0 -> result overflow
;-----------------------------------------------------------------------
proc gg_d32a ; Д32А
;-----------------------------------------------------------------------
mov a,l
sub c
mov a,h
sbb b
rnc
call negb
xra a
.loop:
dad h
rar
xchg
dad h
xchg
jnc .m1
inx h
.m1:
push h
dad b
jnc .m2
ral
.m3:
inx d
inx sp
inx sp
adi 16
jnc .loop
; stc
ret
.m2:
ral
jc .m3
pop h
adi 16
jnc .loop
; stc
ret
;-----------------------------------------------------------------------
endp
;-----------------------------------------------------------------------
proc negb
mov a,b
cma
mov b,a
mov a,c
cma
mov c,a
inx b
ret
endp
;-----------------------------------------------------------------------

[b2m]

Моё изобретение :) Столбиком в двоичном виде:

SHFT: XCHG
DAD H
XCHG
DAD H
RNC
INX D
RET

; вход: HL - квадрат числа
; выход: B - число
; DE - разница между исходным числом и квадратом
SQRT: LXI B,8
LXI D,0
SQRT1: CALL SHFT
CALL SHFT
PUSH H
MOV A,B
ADD A
MOV B,A
CMA
MOV L,A
MVI H,0FFh
DAD H
INX H
DAD D
JNC SQRT2
INR B
XCHG
SQRT2: POP H
DCR C
JNZ SQRT1
RET

Всего 8 итераций!

[shoorick]

работает!
принимаем на вооружение ;)

[shoorick]

А заодно в виде упражнения и кубический 8)

;================================================= ======================
; E = CBRT(HL)
;-----------------------------------------------------------------------
cbrt:
;-----------------------------------------------------------------------
push h
lxi d,16 shl 8
.loop:
mov a,e
add e
mov e,a
mov h,e
inr h
push d
call umul88
pop d
push h
pop b
dad h
dad b
inx h
mov a,d
dcr a
jz .ok
.shift:
dad h
jc .skip
dcr a
jnz .shift
.ok:
pop b
mov a,c
sub l
mov l,a
mov a,b
sbb h
jc .endif
mov b,a
mov c,l
inr e
.endif:
push b
.skip:
mov a,d
sui 3
mov d,a
jp .loop
pop b
ret
;================================================= ======================

использует умножение из той же книжки :)

;-----------------------------------------------------------------------
; HL = H * E
; BС,E:CONST
;-----------------------------------------------------------------------
umul88 equ gg_u88b1
;-----------------------------------------------------------------------
proc gg_u88b1 ; У88Б1 модифицированная
;-----------------------------------------------------------------------
mvi d,0
mvi l,0
mvi a,8
.loop:
dad h
jnc .skip
dad d
.skip:
dcr a
jnz .loop
ret
;-----------------------------------------------------------------------
endp
без него не смог, так как портировал данный код из другой известной книжки:

int icbrt(unsigned x){
int s;
unsigned y,b;
s=15;
y=0;
while (s>=0) {
y = 2*y;
b = (3*y*(y+1)+1)<<s;
//b = b & 0xFFFF;
if (x>=b) {
x = x - b;
y = y + 1;
}
s = s - 3;
//printf("x=%6u, y=%6u, b=%6u, s=%2i\n",x,y,b,s);
}
return y;
}
:)

[b2m]

портировал данный код из другой известной книжки
По структуре очень похож на то, что я сделал выше для квадратного. Наверное, несложно будет модифицировать моё творение, чтобы оно кубический корень извлекало.

Вот моё в таком-же виде:

int isqrt(unsigned x){
int s;
unsigned y,b;
s=14;
y=0;
while (s>=0) {
y = 2*y;
b = (2*y+1)<<s;
if (x>=b) {
x = x - b;
y = y + 1;
}
s = s - 2;
}
return y;
}

- - - Добавлено - - -

Только вот вычислять 3*y*(y+1)+1 вместо 2*y+1 чуточку сложнее :)

- - - Добавлено - - -

Интересно, а можно вывести формулу для b = ???, чтобы оно например корень 5-й степени извлекало?

- - - Добавлено - - -

2*y+1 = (y+1)^2 - y^2
3*y*(y+1)+1 = 3*y^2 + 3*y + 1 = (y+1)^3 - y^3
Что у нас там для пятой степени получается? :)

- - - Добавлено - - -

(y+1)^5 - y^5 = 5*y^4 + 10*y^3 + 10*y^2 + 5*y + 1 = 5*y*(y*(y*(y+2)+2)+1)+1
Щас проверим :)

- - - Добавлено - - -

Ха, работает :)

int isqrt5(unsigned x){
int s;
unsigned y,b;
s=20;
y=0;
while (s>=0) {
y = 2*y;
b = (5*y*(y*(y*(y+2)+2)+1)+1)<<s;
if (x>=b) {
x = x - b;
y = y + 1;
}
s = s - 5;
}
return y;
}

Разрядность для тестов увеличил.

[shoorick]

да, не попалась мне такая sqrt... :)
в хакер-делайт была такая:
int isqrt(unsigned x) {
unsigned m,y,b;
m=0x40000000;
y = 0;
while (m!=0) {
b = y | m;
y = y >> 1;
if (x>=b) {
x = x - b;
y = y | m;
}
m = m >> 2;
}
return y;

}

-- ввела меня немного в ступор :D

а с разрядностью - да, попотел, пока понял, что при значениях > 37^3 переполнялось b, победил при помощи "jc .skip"
надо будет попробовать 32-битные процедуры замутить :) как минимум sqrt

[b2m]

ввела меня немного в ступор
Да, правильность не очевидна :) И что, работает?

- - - Добавлено - - -

Хотя да, практически то-же самое. Только сдвиги в другую сторону.

[shoorick]

да, работает :)

[iceoflame]

Напоминает шутку:

ya_frosia: Связка ломов, как правило, тонет
alexei: ya_frosia: Но в ртути прекрасно плавает
zoogenic: alexei: Но если ломы урановые, то и во ртути тонут
alexei: zoogenic: сам топи урановые ломы в ртути

https://lurkmore.co/%D3%F0%E0%ED%EE%E2%FB%E5_%EB%EE%EC%FB

[shoorick]

в вопросе материала для ломов главное помнить о кейворине!

[shoorick]

32-bit

;================================================= ======================
; BC = SQRT(DEHL)
;-----------------------------------------------------------------------
sqrt32:
;-----------------------------------------------------------------------
mvi a,30 + 2
push psw ; s+2
lxi b, 0 ; y
push d ; xh
push h ; xl
;-----------------------------------------------------------------------
.loop:
pop h
pop d
pop psw
sui 2 ; a=s
rm
;-----------------------------------------------------------------------
push psw
push d
push h
mov l,c
mov h,b
dad h
mov c,l
mov b,h ; BC=y=2*y
;-----------------------------------------------------------------------
lxi d,0
dad h
inx h ; DEHL=b=2*y+1
;-----------------------------------------------------------------------
inr a
.shift:
dcr a
jz .skip
xchg
dad h
jc .loop
xchg
dad h
jnc .shift
inx d
jmp .shift
.skip: ; DEHL << A=s
;-----------------------------------------------------------------------
push d
push h
lxi h,4
dad sp
xchg ; DE --> X
lxi h,0
dad sp ; HL --> B
;-----------------------------------------------------------------------
ldax d
sub m
mov m,a
inx h
inx d
ldax d
sbb m
mov m,a
inx h
inx d
ldax d
sbb m
mov m,a
inx h
inx d
ldax d
sbb m
mov m,a
;-----------------------------------------------------------------------
pop h
pop d
jc .loop
;-----------------------------------------------------------------------
xthl
pop h
xchg
xthl
push d
inx b
jmp .loop
;================================================= ======================

можно и посмотреть: 57758

[ivagor]

Вариант b2ma удачно сочетает компактность и скорость. Если добавить 3 байта, то будет на 15% быстрее. Но если действительно нужно много и быстро (почти в 3 раза быстрее) считать квадратные корни, то пригодится вариант с таблицей (перевел с z80 отсюда (http://ivr.webzone.ru/articles/dop_ar/#m5)).

[ivagor]

Немного сокращенная и ускоренная версия табличного варианта (предыдущий убрал)

[b2m]

Если добавить 3 байта, то будет на 15% быстрее.
Ну ты как всегда, в своём репертуаре :)
Колись уже, где там 3 байта добавлять надо.

[HardWareMan]

Вконце, 3 NOPа после RET. :)

[ivagor]

Ну я действительно в своем репертуаре - всего лишь избавился от процедуры

[HardWareMan]

А действительно.

0001 0000 ; вход: HL - квадрат числа
0002 0000 ; выход: B - число
0003 0000 ; DE - разница между исходным числом и квадратом
0004 0000 01 08 00 SQRT: LXI B,8
0005 0003 11 00 00 LXI D,0
0006 0006 CD 22 00 SQRT1: CALL SHFT
0007 0009 CD 22 00 CALL SHFT
0008 000C E5 PUSH H
0009 000D 78 MOV A,B
0010 000E 87 ADD A
0011 000F 47 MOV B,A
0012 0010 2F CMA
0013 0011 6F MOV L,A
0014 0012 26 FF MVI H,0FFh
0015 0014 29 DAD H
0016 0015 23 INX H
0017 0016 19 DAD D
0018 0017 D2 1C 00 JNC SQRT2
0019 001A 04 INR B
0020 001B EB XCHG
0021 001C E1 SQRT2: POP H
0022 001D 0D DCR C
0023 001E C2 06 00 JNZ SQRT1
0024 0021 C9 RET
0025 0022 ;
0026 0022 EB SHFT: XCHG
0027 0023 29 DAD H
0028 0024 EB XCHG
0029 0025 29 DAD H
0030 0026 D0 RNC
0031 0027 13 INX D
0032 0028 C9 RET
Против:

0001 0000 ; вход: HL - квадрат числа
0002 0000 ; выход: B - число
0003 0000 ; DE - разница между исходным числом и квадратом
0004 0000 01 08 00 SQRT: LXI B,8
0005 0003 11 00 00 LXI D,0
0006 0006 EB SQRT1: XCHG
0007 0007 29 DAD H
0008 0008 EB XCHG
0009 0009 29 DAD H
0010 000A D2 0E 00 JNC SQRT11
0011 000D 13 INX D
0012 000E EB SQRT11: XCHG
0013 000F 29 DAD H
0014 0010 EB XCHG
0015 0011 29 DAD H
0016 0012 D2 16 00 JNC SQRT12
0017 0015 13 INX D
0018 0016 E5 SQRT12: PUSH H
0019 0017 78 MOV A,B
0020 0018 87 ADD A
0021 0019 47 MOV B,A
0022 001A 2F CMA
0023 001B 6F MOV L,A
0024 001C 26 FF MVI H,0FFh
0025 001E 29 DAD H
0026 001F 23 INX H
0027 0020 19 DAD D
0028 0021 D2 26 00 JNC SQRT2
0029 0024 04 INR B
0030 0025 EB XCHG
0031 0026 E1 SQRT2: POP H
0032 0027 0D DCR C
0033 0028 C2 06 00 JNZ SQRT1
0034 002B C9 RET

[shoorick]

т.е. табличный быстрее, просто жрет OVER 512 байт памяти.

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

думаю, оптимальный вариант будет зависить от конкретных условий :) гдавное, есть из чего выбирать ;)

[AndTorp]

Использует подпрограмму деления Д32А из книжки Гуртовцева и Гудыменко

портировал данный код из другой известной книжки
shoorick, из какой "другой известной книжки"?

[shoorick]

Hacker's Delight (по-нашему: Алгоритмические трюки для программистов)

[ivagor]

Стоит ли выкладывать очень быстрый, но приближенный вариант или здесь только 100% точные (с усечением)? Сделал по этому (http://www.zxpress.ru/article.php?id=3607&lng=eng) алгоритму. Приближенность связана с тем, что используются только 8 бит мантиссы исходного числа. Но погрешность результата там максимум 1.
Самый быстрый вариант требует аж 7424 байта таблиц, зато (без call и ret) выполняется за 74 такта. На z80 за 69 тактов (!!!), т.е. быстрее, чем авторский вариант. Кроме того, пока разбирался, сделал несколько вариантов, в т.ч. с таблицами всего на 768 байт - все равно быстрее sqrttab2.
Можно попытаться приделать корректор, дожимающий точность до 100%

[shoorick]

ну конечно же сто́ит! это же форум для всех - кому-нибудь пригодится, например, звук генерировать или для фильтра цифрового, или графика - там скорость важнее точности, тем более, что отбрасывание дробной части уже не есть точность ;)

[b2m]

Ну я действительно в своем репертуаре - всего лишь избавился от процедуры
Ну вот, половину кода переписал :) А говорил три байта добавил...

[ivagor]

Приделал к "приблизительной" процедуре корректор - теперь дает на 100% точные результаты. В худшем случае 204 такта (это вместе с call и ret) и пришлось добавить таблицу на 512 байт из sqrttab. Для сравнения sqrttab2 в лучшем случае выполняется за 527 тактов. В среднем sqrtfast7 еще быстрее (в sqrtab2 медленнее), разница почти в 4 раза.
Для подавляющего большинства применений можно поставить ret после ldax b и пользоваться приближенным результатом (только обратите внимание, что результат в таком варианте будет в A, не в L), тогда параметры будут как я написал выше.
Само вычисление корня по этому алгоритму вряд ли получится ускорить, а вот корректор скорее всего можно оптимизировать.

- - - Добавлено - - -

Пара слов насчет точности.
В приложенном Etalon.zip две таблицы результатов, рассчитанные в матлабе
sqrtmatlab.fix - с округлением к нулю, как считают практически все точные целочисленные извлекалки (для примера проверил sqrttab2 - совпадает).
sqrtmatlab.round - с округлением к ближайшему значению.

Что дает sqrtfast7
Выложенный ранее вариант с корректором - совпадает с sqrtmatlab.fix
Если из sqrtfast7 убрать корректор - получится sqrtfast7.round (из sqrtfast7results.zip)
Приложенный к данному посту sqrtfast7fix.zip - выдает sqrtfast7.fix (тоже из sqrtfast7results.zip)
Чтобы оценить точность приближенных вариантов можно сравнить sqrtfast7.round с sqrtmatlab.round, а sqrtfast7.fix с sqrtmatlab.fix.
Кстати, корректор от предыдущего варианта к sqrtfast7fix не совсем подходит, прицеплять его туда смысла нет.

[ivagor]

На базе варианта b2mа и вот этого (http://www.zxpress.ru/zxnet/code.zx/658) варианта для z80 сделал оптимизированную версию без таблиц.
sqrtnotab - в 3 раза больше варианта b2mа (121 байт) и в 3.0-3.2 раза быстрее
sqrtnotabCompact - в 2 раза больше варианта b2mа (83 байта) и в 2.6-2.8 раза быстрее

sqrtnotab догнал sqrttab2, который с его таблицей 512 байт теперь точно не нужен. sqrtfast7 конечно все равно быстрее раза в 3-3.5 (это если считать точно, а если приближенно то в 6 раз), но и памяти жрет немерено.
Возможно вычисление двух старших бит еще можно оптимизировать.

[shoorick]

вариантов уже больше чем сфер применения :D
не успеваю тестировать 8)
кстати, встречал вариант на Z80 с округлением ответа (не обрезанием),
завтра с работы ссылку покажу.

[shoorick]

вот оно: http://z80-heaven.wikidot.com/math#toc27
но, к сожалению, аргумент 8-битный

а это: http://z80-heaven.wikidot.com/math#toc30
вариант с волшебным 400Н, сдвигающимся вправо ;)

[ivagor]

Удалил свой первый вариант

- - - Добавлено - - -


вот оно: http://z80-heaven.wikidot.com/math#toc27
но, к сожалению, аргумент 8-битный
Это ерунда, а не округление. Достаточно одного неверного результата, чтобы нельзя было сказать, что процедура округляет к ближайшему целому, вот для примера парочка (там еще есть):
sqrt(2)=1,41~1; RoundSqrtE(2)=2
sqrt(6)=2,45~2; RoundSqrtE(6)=3

- - - Добавлено - - -

Однако RoundSqrtE можно допилить напильником даже не увеличив размер - меняем в конце
jr c,$+3
на два условных retа
ret c
ret z
и результат совпадет с первыми 256 байтами таблицы sqrtmatlab.round из Etalon.zip

"Сделал как у них, только правильно". Добавил округление к sqrtnotab (первые варианты с округлением по таблице удалил - у них нет преимуществ).
В качестве примера два варианта:
57933 - округление с насыщением (чтобы результат поместился в байт). Результат совпадает с таблицей sqrtmatlab.round, которую я выкладывал (http://zx-pk.ru/threads/26834-kvadratnyj-koren-na-i8080.html?p=882151&viewfull=1#post882151) ранее в Etalon.zip
57932 - "полное" округление без насыщения. Пришлось выделить под результат регистровую пару, хотя из старшего регистра нужен только один бит (да и то редко).

[shoorick]

да, всё надо проверять, даже очевидное...
но, честно говоря, приветствуя скорость, меня смущает разрядность... причем 32-битный и даже 64-битный (33...34...40-битный) аргумент вполне реален и мне не кажется роскошью: например, если в Специалисте провести диагональ на экране, то ее длина будет как раз sqrt(256^2+384^2)=sqrt(212992)=sqrt(34000h)~461.5= 461
хотя в таких случаях можно аргумент "делить" сдвигом на 4 или 16, а затем "умножать" на 2 или 4 соответственно: sqrt(128^2+192^2)=sqrt(53248)=sqrt(D000h)~230.7=23 0*2=460

[ivagor]

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

- - - Добавлено - - -

32 битный вариант. В 4,5 раза больше ньютона (http://zx-pk.ru/threads/26834-kvadratnyj-koren-na-i8080.html?p=880965&viewfull=1#post880965), но раз в 10 быстрее (и еще точно можно оптимизировать). Отмечу, что используется самомодифицирующийся код (этот фрагмент хорошо бы переписать).

[shoorick]

Тем временем реализовал деление на 10 сдвигами:

unsigned divu10(unsigned n) {
unsigned q, r;
q = (n >> 1) + (n >> 2);
q = q + (q >> 4);
q = q + (q >> 8);
q = q + (q >> 16);
q = q >> 3;
r = n - q*10;
return q + ((r + 6) >> 4);
// return q + (r > 9);
} -- отсюда (http://www.hackersdelight.org/divcMore.pdf)


;================================================= ======================
; HL [A] = HL/10
;-----------------------------------------------------------------------
proc udiv10
;-----------------------------------------------------------------------
ora a
mov a,h
rar
mov b,a
mov a,l
rar
mov c,a ; bc = n >> 1
ora a
mov a,b
rar
mov d,a
mov a,c
rar
add c
mov e,a
mov c,a
mov a,d
adc b
mov b,a
mov d,a ; bc = de = (hl>>1)+(hl>>2) = q
xchg ; de = n ; bc = hl = q
xra a
call shr_hl.s4 ; hl = q >> 4
dad b ; hl = q = q + (q>>4)
mov c,h
mvi b,0
dad b ; hl = q = q + (q>>8)
xra a
call shr_hl.s3 ; hl = q = q >> 3
mov b,h
mov c,l ; bc = q
dad h
dad h
dad b
dad h ; hl = q*10
mov a,e
sub l
adi 6
rar
rar
rar
rar
ani 15 ; a = (r+6)>>4
add c
mov l,a
mov a,b
aci 0
mov h,a ; hl = q + ((r+6)>>4) = n/10
mov a,l
add a
add a
add l
add a ; a = l*10
sub e
cma
inr a ; a = e - l*10 (remainder)
ret
;-----------------------------------------------------------------------
.dummy = shr_hl
;-----------------------------------------------------------------------
endp
;================================================= ======================
;
;================================================= ======================
; FOR .s1-.s7 A MUST BE 0 !
;-----------------------------------------------------------------------
proc shr_hl
;-----------------------------------------------------------------------
xra a
.s1:
dad h
ral
.s2:
dad h
ral
.s3:
dad h
ral
.s4:
dad h
ral
.s5:
dad h
ral
.s6:
dad h
ral
.s7:
dad h
ral
mov l,h
mov h,a
ret
;-----------------------------------------------------------------------
endp
;================================================= ======================
-- надо будет ускорить вывод чисел
++++++++++++++++++++++++++
оказалось, что исходный вариант делит хорошо, но остаток у него левый ;) причем не всегда, а когда аргумент равен 10*х и 10*х+1
поэтому доточил в конце вычисление правильного остатка ;) уж даже не знаю, но надеюсь, что быстрее, чем просто деление на 10
(уже влом такты считать)

[shoorick]

32-битный кубический корень :)
как обычно взял исходник на си из тех же мест, не смог обойтись стеком - использовал память,
доперепилил ;)

int icbrt2(unsigned x) {
int s;
unsigned y, b, bs, y2;

y2 = 0;
y = 0;
for (s = 30; s >= 0; s = s - 3) {
y2 = 4*y2;
y = 2*y;
b = 3*(y2 + y) + 1;
bs = b << s;
if (x >= bs && b == (bs >> s)) {
x = x - bs;
y2 = y2 + 2*y + 1;
y = y + 1;
}
}
return y;
}

это процедурка
;================================================= ======================
; BC = CBRT(DEHL)
;-----------------------------------------------------------------------
cbrt32:
;-----------------------------------------------------------------------
shld .x
xchg
shld .x+2
lxi b,0
mov h,b
mov l,c
shld .y2
shld .y2+2
mvi a,30
;-----------------------------------------------------------------------
.loop:
sta .s
;-----------------------------------------------------------------------
lhld .y2+2
xchg
lhld .y2
mvi a,2
call shl_dehl_a
shld .y2
xchg
shld .y2+2 ; y2=4*y2
xchg
mov a,c
add a
mov c,a
mov a,b
adc a
mov b,a ; y=2*y
dad b
jnc .m1
inx d
.m1: ; dehl = y2 + y
xra a
push b
mov b,h
mov c,l
dad b
aci 0
dad b
aci 0
xchg
mov b,h
mov c,l
dad b
dad b
add l
mov l,a
mov a,h
aci 0
mov h,a
xchg
inx h
mov a,l
ora h
jnz .m2
inx d
.m2:
pop b ; dehl = b = 3*(y2 + y)+1
lda .s
call shl_dehl_a ; dehl = bs = b << s
lda .x
sub l
mov l,a
lda .x+1
sbb h
mov h,a
lda .x+2
sbb e
mov e,a
lda .x+3
sbb d
mov d,a ; dehl = x - bs
jc .skip
shld .x ; x = x - bs
xchg
shld .x+2
lhld .y2+2
xchg
lhld .y2
dad b
jnc .m3
inx d
.m3:
dad b
jnc .m4
inx d
.m4:
inx h
mov a,l
ora h
jnz .m5
inx d
.m5:
shld .y2
xchg
shld .y2+2 ; y2 = y2 + 2*y + 1
inx b ; y = y + 1
.skip:
lda .s
sui 3
jnc .loop
;-----------------------------------------------------------------------
ret
;-----------------------------------------------------------------------
.x dd ?
.y2 dd ?
.s db ?
;================================================= ======================

а это - служебная сдвигайка влево
;================================================= ======================
; DEHL = DEHL << A
;-----------------------------------------------------------------------
proc shl_dehl_a
;-----------------------------------------------------------------------
ora a
rz
;-----------------------------------------------------------------------
cpi 32
jnc .zero_dehl
push psw
ani 7
inr a
.loop:
dcr a
jz .move
xchg
dad h
xchg
dad h
jnc .loop
inx d
jmp .loop
.move:
pop psw
ani 18h
rar
rar
rar
rar
jnc .next
mov d,e
mov e,h
mov h,l
mvi l,0
.next:
rar
rnc
xchg
lxi h,0
ret
.zero_dehl:
lxi d,0
lxi h,0
ret
;-----------------------------------------------------------------------
endp
;================================================= ======================

Приблизительный 64-битный квадратный корень.
Использует указанные выше процедуры.
При необходимости сначала побайтно сдвигает аргумент, а потом соответственно и результат, поэтому для аргументов >32 бит является приблизительным.

;================================================= ======================
; DEHL = SQRT([HL]) APPROX FOR > 32 BIT (VIA 32-BIT SQRT)
;-----------------------------------------------------------------------
proc sqrt64a
;-----------------------------------------------------------------------
lxi b,7
dad b
mvi b,4
xra a
.loop:
cmp m
jnz .ok
dcx h
dcr b
jnz .loop
.ok:
push b
mov d,m
dcx h
mov e,m
dcx h
mov a,m
dcx h
mov l,m
mov h,a
call sqrt32
lxi d,0
mov h,b
mov l,c
pop psw
add a
rz
add a
jmp shl_dehl_a
;-----------------------------------------------------------------------
endp
;================================================= ======================