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

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

Код:

---

;-----------------------------------------------------------------------
;  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   
;-----------------------------------------------------------------------

---

[свернуть]

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

Код:

---

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 итераций!

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

А заодно в виде упражнения и кубический 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;
}

---

[свернуть]

:)

Цитата:

---

Сообщение от shoorick

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

---

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

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

Код:

---

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;
}

---

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

да, не попалась мне такая 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

Цитата:

---

Сообщение от shoorick

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

---

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

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

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

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

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

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

https://lurkmore.co/%D3%F0%E0%ED%EE%...5_%EB%EE%EC%FB

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

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

Вариант b2ma удачно сочетает компактность и скорость. Если добавить 3 байта, то будет на 15% быстрее. Но если действительно нужно много и быстро (почти в 3 раза быстрее) считать квадратные корни, то пригодится вариант с таблицей (перевел с z80 отсюда).

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

Цитата:

---

Сообщение от ivagor

Если добавить 3 байта, то будет на 15% быстрее.

---

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

Вконце, 3 NOPа после 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 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

---

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

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

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

Цитата:

---

Сообщение от shoorick

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

---

Цитата:

---

Сообщение от shoorick

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

---

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

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

Стоит ли выкладывать очень быстрый, но приближенный вариант или здесь только 100% точные (с усечением)? Сделал по этому алгоритму. Приближенность связана с тем, что используются только 8 бит мантиссы исходного числа. Но погрешность результата там максимум 1.
Самый быстрый вариант требует аж 7424 байта таблиц, зато (без call и ret) выполняется за 74 такта. На z80 за 69 тактов (!!!), т.е. быстрее, чем авторский вариант. Кроме того, пока разбирался, сделал несколько вариантов, в т.ч. с таблицами всего на 768 байт - все равно быстрее sqrttab2.
Можно попытаться приделать корректор, дожимающий точность до 100%

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

Цитата:

---

Сообщение от 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 не совсем подходит, прицеплять его туда смысла нет.

На базе варианта b2mа и вот этого варианта для 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 раз), но и памяти жрет немерено.
Возможно вычисление двух старших бит еще можно оптимизировать.

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

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

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

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

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

Цитата:

---

Сообщение от shoorick

вот оно: 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, которую я выкладывал ранее в Etalon.zip
Вложение 57932 - "полное" округление без насыщения. Пришлось выделить под результат регистровую пару, хотя из старшего регистра нужен только один бит (да и то редко).

да, всё надо проверять, даже очевидное...
но, честно говоря, приветствуя скорость, меня смущает разрядность... причем 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

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

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

32 битный вариант. В 4,5 раза больше ньютона, но раз в 10 быстрее (и еще точно можно оптимизировать). Отмечу, что используется самомодифицирующийся код (этот фрагмент хорошо бы переписать).

Тем временем реализовал деление на 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);
}

---

-- отсюда

Скрытый текст

Код:

---

;=======================================================================
;  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
(уже влом такты считать)

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
;=======================================================================

---

[свернуть]