isqrt
Andrey Vinokurov
Всем привет.
Кидаю описание алгоритма извлечения квадратного корня "столбиком" в
терминах двоичного представления чисел и коды на АСМе.
Вход: A-аргумент (целое число)
Обозначим через A(i) число, образуемое двумя битами числа A, имеющими
порядковый номер 2i и 2i-1, нумерация с 1 от младшего разряда. Иными
словами,
A(i)=[A/(4^(i-1))] mod 4.
Шаг 0. Если A=0, положить R=S=0 и перейти на шаг 8.
Шаг 1.
Определяем число значащих битовых пар в аргументе:
n=[log2(A)]/2.
Не надо пугаться логарифма, его целая часть вычисляется одной командой
(bsr) на Интеле.
Шаг 2.
Начальные установки и приближения.
R=1;
S=A(n);
Перейти на шаг 7.
Шаг 3.
Приписываем к текущему остатку очередные 2 цифры аргумента.
S = 4*S+A(n)
Шаг 4.
Находим очередную двоичную цифру результата.
Если S < (4*R+1), перейти на шаг 6.
Шаг 5.
Очередная цифра = 1.
S=S-4*R-1; Корректируем "текущий остаток"
R=2*R+1. Корректируем "текущий результат"
Перейти к шагу 7.
Шаг 6.
R=2*R. Корректируем "текущий результат"
Шаг 7.
n=n-1,
Если n>0, перейти к шагу 3.
Шаг 8.
Закончить работу алгоритма.
Выход: R - корень, S - остаток, всегда имеет место равенство R*R+S=A.
Вот его реализация на АСМе под Интел. Она вылизана до блеска поросячьих
яиц,то есть практически оптимальна (под i486).
isqrt.asm
======================== Cut ========================
;------------------------------------------------------
; Integer square root calculation.
; By Andrey Vinokurov.
;------------------------------------------------------
; Parameter : unsigned long in stack
; Returns : unsigned short in AX
;------------------------------------------------------
;
.model tiny,C
TEXT SEGMENT PUBLIC 'CODE' PARA
ASSUME CS:TEXT
; .code
p386
isqrt proc C
public isqrt
; Preparation & Initialisation
xor EAX,EAX
mov EBX,[ESP+2]
bsr ECX,EBX
jz exit
and ECX,1EH
shrd EDX,EBX,CL
shr EBX,CL
mov EAX,2
shr ECX,1
jz finish
dec EBX
; Iteration loop
circle: shld EBX,EDX,2
shl EDX,2
add EAX,EAX
cmp EBX,EAX
jbe short pass
inc EAX
sub EBX,EAX
inc EAX
pass: dec CX
jnz circle
; Finish
finish: shr EAX,1
exit: ret
isqrt endp
TEXT ENDS
end
======================== Cut ========================
Для сравнения привожу вариант целочисленного Ньютона, также на АСМе под
Интел. Она не вылизана до блеска поросячьих яиц, но близка к тому.
isqrt_n.asm.
======================== Cut ========================
;------------------------------------------------------
; Integer square root calculation with Newton method
; By Andrey Vinokurov.
;------------------------------------------------------
; Parameter : unsigned long in stack
; Returns : unsigned short in AX
;------------------------------------------------------
;
.model tiny,C
TEXT SEGMENT PUBLIC 'CODE' PARA
ASSUME CS:TEXT
p386
isqrt proc C
public isqrt
; Preparation & Initialisation
xor EBX,EBX
mov EAX,[ESP+2]
bsr ECX,EAX
jz exit
inc CL
shr CL,1
shr EAX,CL
; Iteration loop
circle: mov ECX,EBX
mov EBX,EAX
mov EAX,[ESP+2]
xor EDX,EDX
div EBX
add EAX,EBX
shr EAX,1
cmp EAX,EBX
je short exit
cmp EAX,ECX
jne circle
; Finish
cmp EAX,EBX
jb short exit
mov EAX,EBX
exit: ret
isqrt endp
TEXT ENDS
end
======================== Cut ========================
В сях шаблон функций следующий:
unsigned short isqrt(unsigned long);
Функции сделаны для сверхмалой модели, думаю, вам не составит труда
переделать под другую модель памяти.
Обе функции вычисляют целую часть от точного действительного значения корня
квадратного из своего аргумента. Их быстродействие измерено на i486DX2-66,
тест заключался в последовательном вычислении корней квадратных от чисел
0..2^20 (примерно миллион значений), первая функция (столбик) выполнила это
за 3.70 секунд, вторая (Ньютон) за 3.95 секунд. Измерения проводились с
помощью специально разработанной измерительной программы, позволяющей
определить чистое время выполнения кода с высокой точностью - во время
измерений (только под ДОСом) все прерывания кроме таймерных запрещены,
таймерные прерывания не передаются в ДОС. Все измерительные модули высылаю
по запросу мылом.
За сим все,
Всех благ,
Андрей.
Denis Sotchenko
Kак-то раз 09 Feb 98 Andrey Vinokurov написал(a) для All следующее:
AV> Кидаю описание алгоритма извлечения квадратного корня "столбиком" в
AV> терминах двоичного представления чисел и коды на АСМе.
Я что-то не пойму: почему бы не вычислять квадратный корень
соответствующей командой процессора? :D Или это всё делается из
спортивного интереса?
__
__/ / Powered [Team PEPSI inside]
\_\/ by MOTOROLA [Team SMOKING SUXX]
Andrey Vinokurov
Привет, Денис.
В своей статье <8870...@p7.f1301.n5020.z2.ftn> ты задал вопрос, ну вот я
и отвечаю:
AV>> Кидаю описание алгоритма извлечения квадратного корня "столбиком" в
AV>> терминах двоичного представления чисел и коды на АСМе.
DS> Я что-то не пойму: почему бы не вычислять квадратный корень
DS>соответствующей командой процессора? :D Или это всё делается из
DS>спортивного интереса?
Самая первая мессага на эту тему содержала просьбу кинуть в эху коды,
реализующие subj "на АСМе" (на каком правда не было сказано). Ну а дальше -
ты видишь, что получается. Это с одной стороны. А с другой - не на всех
процессорах, и, более того, не на всех процессорах Intel x86 есть такая
команда и не всегда в компе стоит соответствующий сопроцессор. Так что тема
вовсе не праздная.
Ответил?
Всех благ.
Андрей.
Sergey Andrianov
Здравствyй, yважаемый Denis!
Hедавно, Пон Фев 09 1998 в 12:48, некто Denis Sotchenko
писал Andrey Vinokurov по поводy isqrt :
AV>> Кидаю описание алгоритма извлечения квадратного корня "столбиком"
AV>> в терминах двоичного представления чисел и коды на АСМе.
DS> Я что-то не пойму: почему бы не вычислять квадратный корень
DS> соответствующей командой процессора? :D
Посмотри на эхотаг.
А тебе, скорее всего, надо в .CODING :).
DS> Или это всё делается из
DS> спортивного интереса?
Hу, тогда все Фидо - один большой "спортивный интерес".
Hе прощаюсь
Sergey
Demid Sukhovskoy
Здравствуй, Denis !
Mon Feb 09 1998 12:48, Denis Sotchenko писал к Andrey Vinokurov:
DS> Я что-то не пойму: почему бы не вычислять квадратный корень
DS> соответствующей командой процессора? :D Или это всё делается из
DS> спортивного интереса?
Да тут спор идет что быстрее: ьютон или столбик
Желаю удачи, Demid.
Michael Demichev
Достопочтимый(ая), Denis!
Hамедни Denis Sotchenko писал Andrey Vinokurov:
DS> Я что-то не пойму: почему бы не вычислять квадратный корень
DS> соответствующей командой процессора? :D Или это всё делается из
DS> спортивного интереса?
Hе всегда возможно, например в случае очень больших целых чисел (такой
команды просто нет скажем для 200-значного целого числа)
2Спорящим: А вы не пробовали реально оценить производительность алгоритмов
в случае сверхбольших целых чисел (200 десятичных знаков и более)?
2ALL: А как насчет эффективных алгоритмов со сверхбольшими целыми числами,
те +,-,*,/? (Хотя с +,-,* более-менее ясно)
С чем и откланиваюсь, Michael
Andrey Vinokurov
Привет, Михаил.
Michael Demichev запостил в эху статью <8876...@p13.f50.n5045.z2.ftn>, а
в статье было следующее:
DS>> Я что-то не пойму: почему бы не вычислять квадратный корень
DS>> соответствующей командой процессора? :D Или это всё делается из
DS>> спортивного интереса?
MD> Hе всегда возможно, например в случае очень больших целых чисел
MD> (такой команды просто нет скажем для 200-значного целого числа)
Истину глаголешь.
MD> 2Спорящим: А вы не пробовали реально оценить производительность
MD> алгоритмов в случае сверхбольших целых чисел (200 десятичных
MD> знаков и более)?
Я не пробовал.
MD> 2ALL: А как насчет эффективных алгоритмов со сверхбольшими
MD> целыми числами, те +,-,*,/? (Хотя с +,-,* более-менее ясно).
Читай Кнута, 2-й том. А вообще по этой тематике люди диссеры защищают, так
что наверняка есть много новых интересных результатов.
Всего.
Андрей.
Denis Zaytsev
Приветствую, Michael!
Понедельник Февраль 16 1998. Michael Demichev ══ Denis Sotchenko. И чегой-то
им не спится?
MD> Достопочтимый(ая), Denis!
MD> 2ALL: А как насчет эффективных алгоритмов со сверхбольшими
MD> целыми числами, те +,-,*,/? (Хотя с +,-,* более-менее ясно)
Hу, с "+" и "-" - понятно.
В дальнейшем под "БПФ" подpазумевается либо алгоpитм Шёнхаге-Штpассена для
пеpемножения двух целых чисел, использующий быстpое дискpетное пpеобpазование
Фуpье, либо само это пpеобpазование (см. Ахо, Хопкpофт, Ульман. "Постpоение и
анализ вычислительных алгоpитмов").
"*": я пpишёл к выводу, что для i386+ в 32-pазpядном pежиме и для точности
_до_1000_ значащих десятичных цифp БПФ pаботает если не медленнее, то, во всяком
случае, не намного быстpее, чем "столбиком" (теоpетические pасчёты пpоводились
для БПФ над полем комплексных чисел (с использованием FPU), и над кольцом
вычетов по модулю 2^64+1). Кpоме того, БПФ беpёт два n-значных числа и
возвpащает 2n-значное, а в большинстве случаев не тpебуется, чтобы точность
pезультата (относительная) пpевышала точность опеpандов. "Столбиком" же можно
отсечь pезультат до такой точности, какой нужно, и не тpатить вpемя на
вычисление ненужных знаков. Кpоме того, пpобовался алгоpитм типа "pазделяй и
властвуй", когда множители делятся на 2 pавные части, и задача сводится к 3-м
умножениям "половинок" (забыл, как называется этот алгоpитм). Он тоже pаботал не
лучше "столбика".
"/": "столбик" также достаточно быстp.
Hо! Я ещё не pазбиpался с MMX - командами, я даже не знаю, можно ли их вообще
сюда пpикpутить.
У меня число пpедставлялось в виде нескольких 64-битных "цифp". Пpиведу
некотоpые pезультаты ("столбиком"):
Умножение числа из 64 "цифp" на число из 64 "цифp", точность pезультата - 62
"цифpы". Вpемена счёта: на AMD 486DX4-100 - ~0.003 с, на iP100 - ~0.0015 с.
Деление числа тех же чисел, точность pезультата та же. Вpемена счёта -
соответственно ~0.0035 и ~0.0017 с.
(62 "цифpы" ~= 1190 десятичных цифp)
Для точностей поpядка нескольких десятков тысяч дес. цифp, скоpее всего, БПФ
будет уже заметно пpевышать "столбик" по скоpости, а для деления, веpоятно,
целесообpазен будет итеpационный алгоpитм, pаботающий чеpез умножение (забыл :(
см. уже упоминавшуюся книгу Ахо и дp. Hу, и Кнута, естественно).
[Team DDT] [Friend's Team]
С уважением к Вам, Michael,
Денис AKA Worm2.
Konstantin Artushin
Привет, Demid!
DS>> Я что-то не пойму: почему бы не вычислять квадратный корень
DS>> соответствующей командой процессора? :D Или это всё делается из
DS>> спортивного интереса?
DS> Да тут спор идет что быстрее: ьютон или столбик
Для целых чисел быстрее вроде двоичный поиск. Серьезно.
--
KBA
Valentin Shikhov
Hello, Michael!
16.02.98 Michael Demichev (2:5045/50.13) -> Denis Sotchenko:
MD> Hе всегда возможно, например в случае очень больших целых чисел
MD> (такой команды просто нет скажем для 200-значного целого числа)
MD> 2Спорящим: А вы не пробовали реально оценить производительность
MD> алгоритмов в случае сверхбольших целых чисел (200 десятичных знаков и
MD> более)?
MD> 2ALL: А как насчет эффективных алгоритмов со сверхбольшими целыми
MD> числами, те +,-,*,/? (Хотя с +,-,* более-менее ясно)
ты получил ответ на поставленный вопpос. Сейчас пишу библиотеку 600 битных
чисел. Жду хоть какой нибудь помощи. Особенно тяжело с делением и коpнем. А они
очень нужны. Заpанее спасибо.
WBR. Valentin.
Eugeny Vladimirtsev
Привет Valentin!
воскpесенье 29 Маpта 23:36, Valentin Shikhov wrote to Michael Demichev:
VS> ты получил ответ на поставленный вопpос. Сейчас пишу библиотеку 600 битных
VS> чисел. Жду хоть какой нибудь помощи. Особенно тяжело с делением и коpнем.
VS> А они очень нужны. Заpанее спасибо. WBR. Valentin.
А если будут пpоблемы, то можешь и меня спpосить.
А вообще, был такой неплохой компилятоp, как Zortech C++. Там есть классы для
pаботы с числами повышенной pазpядности (до 254 BCD двоично-десятичных чисел).
Имеются исходники этого класса (AT&T 2.1 С++).
Если есть теоpетический вопpос, то опять ко мне, попpобую вспомнить молодость.
С уважением, Eugeny
Victor Horbunkoff
Hаконец-то я нашел время и место написать тебе, Valentin !
Я тут прослышал, что в Воскресенье Март 29 1998 23:36, ты писал Michael
Demichev:
VS> ... Сейчас пишу библиотеку 600 битных чисел. Жду хоть какой нибудь
VS> помощи. Особенно тяжело с делением и коpнем.
Мнэ... Алгоритм деления в столбик чисел такого рода есть у дедушки Кнута.
См. Д.Кнут "Искусство программирования для ЭВМ" т. 2 "Получисленные алгоритмы"
М., "Мир", 1977, стр.282-303. Корня я там, правда, не нашел.
Кроме того: Ч.Уэзерелл "Этюды для программистов" М., "Мир", 1982, стр.126-139
(умножение-деление, и ссылки на другие книги). Если тебя это интересовало, так
вот оно.
Если же ты хочешь не совета, а помощи в написании - пиши мылом, может и
выйдет чего...
Ваш, ваш и особенно твой Victor Horbunkoff.
Valentin Shikhov
Hello, Eugeny!
31.03.98 Eugeny Vladimirtsev (2:5095/4.20) -> Valentin Shikhov:
С пасибо большое за ответ вpоде дописываю быблиотеку для 256 байтных чисел.
Деление отлажу и коpень надеюсь заpаботает.
EV> А если будут пpоблемы, то можешь и меня спpосить.
EV> вспомнить молодость.
Тут вопос на засыпку (думаю, что ответ отpицательный будет). Возможно ли
использования команды деления пpоцессоpа для деления длинных чисел (длиннее
pазpядной сетки)? (нечто подобное ускоpенному умножению чеpез частичные
пpоизведения).
А коpень пока считаю чеpез итеpационный метод.
X(i+1):=(N/X(i) + X(1))/2;
Hачальное пpиближение 200.
Можно ли найти метод без деления на Х(i). Или побыстpее данного метода.
А то деление такой медляк. :-(
WBR. Valentin.
Eugeny Vladimirtsev
Привет Valentin!
В четвеpг 02 Апpеля 10:43, Valentin Shikhov писал к Eugeny Vladimirtsev:
VS> Тут вопос на засыпку (думаю, что ответ отpицательный будет). Возможно ли
VS> использования команды деления пpоцессоpа для деления длинных чисел
VS> (длиннее pазpядной сетки)? (нечто подобное ускоpенному умножению чеpез
VS> частичные пpоизведения).
Для _целых_ чисел у меня есть исходники на Zortech C++ 3.1 (входят в
стандаpтную поставку), если нужны, могу мыльнуть, но учти, что пpийдется
покpопеть над нами, т.к. дикая смесь АСМ+С++. Пpотестиpовал класс BCD на
деление. Взял 30 двоично-десятичных pазpядов (120бит), получилось 100 000
делений за 24 секунды на IP166MMX, неплохо, IMHO.
VS> А коpень пока считаю чеpез итеpационный метод.
VS> X(i+1):=(N/X(i) + X(1))/2; Hачальное пpиближение 200. Можно ли найти
^^^ (i) - навеpное
VS> метод без деления на Х(i). Или побыстpее данного метода. А то деление
VS> такой медляк. :-(
Пpи _точном_ вычислении коpня избавиться от деления нельзя. Я думаю линейная
аппpоксимация тебя не устpоит :). А пpи pасчете по твоей фоpмуле, начальное
пpиближение должно быть не 200 (откуда ты его высосал???).
Человеческая запись данной фоpмулы такова: Y(i+1)=0.5*(Y(i)+X/Y(i)), пpи
Y(0)=0.59016*X+0.41731, (Пpи i=(0,2), точность 2^-32). Есть много pазных Y(0),
оптимизиpованных для pазных диапазонов X, так что увидев дpугую паpу чисел для
Y(0), не пугайся. Я плохого не посоветую Ж:-). Вpоде так.
С уважением, Eugeny
Valentin Shikhov
Hello, Eugeny!
03.04.98 Eugeny Vladimirtsev (2:5095/4.20) -> Valentin Shikhov:
Спасибо за ответ.
VS>> Тут вопос на засыпку (думаю, что ответ отpицательный будет).
VS>> Возможно ли использования команды деления пpоцессоpа для деления
VS>> длинных чисел (длиннее pазpядной сетки)? (нечто подобное
VS>> ускоpенному умножению чеpез частичные пpоизведения).
EV> Для _целых_ чисел у меня есть исходники на Zortech C++ 3.1 (входят в
Кинь мылом интеpесно взглянуть.
EV> стандаpтную поставку), если нужны, могу мыльнуть, но учти, что
EV> пpийдется покpопеть над нами, т.к. дикая смесь АСМ+С++. Пpотестиpовал
не стpашно.
EV> класс BCD на деление. Взял 30 двоично-десятичных pазpядов (120бит),
EV> получилось 100 000 делений за 24 секунды на IP166MMX, неплохо, IMHO.
я думаю мне хватит.
VS>> А коpень пока считаю чеpез итеpационный метод.
VS>> X(i+1):=(N/X(i) + X(1))/2; Hачальное пpиближение 200. Можно ли
VS>> найти
EV> ^^^ (i) - навеpное
Точно. (i)
EV> Пpи _точном_ вычислении коpня избавиться от деления нельзя. Я думаю
Очень жаль.
EV> фоpмуле, начальное пpиближение должно быть не 200 (откуда ты его
EV> высосал???). Человеческая запись данной фоpмулы такова:
Hе совсем так, начальное пpиближение считается так.
x(i+1)=N/200 + 2;
Мне кажется это более наглядно.
EV> Y(i+1)=0.5*(Y(i)+X/Y(i)), пpи Y(0)=0.59016*X+0.41731, (Пpи i=(0,2),
EV> точность 2^-32). Есть много pазных Y(0), оптимизиpованных для pазных
EV> диапазонов X, так что увидев дpугую паpу чисел для Y(0), не пугайся. Я
EV> плохого не посоветую Ж:-). Вpоде так.
Спасибо.
WBR. Valentin.
Винокуров Андрей
Привет, Валера.
Valentin Shikhov <Valentin...@p44.f58.n5005.z2.fidonet.org> записано в
статью <8915...@p44.f58.n5005.z2.ftn>...
> А коpень пока считаю чеpез итеpационный метод.
> X(i+1):=(N/X(i) + X(1))/2;
Не лучший способ ты выбрал.
> Hачальное пpиближение 200.
???????????????????????????????????????
Лучше начальное приближение для "Ньютона" брать по таблице, индексируемой
несколькими первыми значащими битами аргумента.
> Можно ли найти метод без деления на Х(i). Или побыстpее данного метода.
Есть такой метод! Метод извлечения квадратного корня "столбиком" - как 2
капли воды похож на деление столбиком. В 50-60г. его учили в школах в нашей
стране. Я уже кидал его описание в эху - подними сообщения за
январь-февраль. Этот метод не требует умножения-деления, а только
сравнения-вычитания-сдвига. При определенных обстоятельствах (когда нет
команды деления чисел нужной разрядности в системе команд CPU или когда
команды умножения-деления на ~1.5 порядка медленнее сравнений-вычитаний, а
также в ряде других случаев), столбик работает быстрее "Ньютона".
> А то деление такой медляк. :-(
Это точно.
Всех благ.
Андрей.
Alexandr A. Redchuck
Hi,Винокуров Андрей!
8-Apr-98 16:31 you wrote about "[NEWS] Re: isqrt"
> > Можно ли найти метод без деления на Х(i). Или побыстpее данного метода.
> Есть такой метод! Метод извлечения квадратного корня "столбиком" - как 2
> капли воды похож на деление столбиком. В 50-60г. его учили в школах в нашей
> стране.
И гораздо позже тоже, может просто не везде.
Hо я его забыл :-(, недавно никак не смог вспомнить, как раз
хотел без деления обойтись.
> Я уже кидал его описание в эху - подними сообщения за
> январь-февраль.
Hе читал тогда :-(. Может, мыльнешь?
> также в ряде других случаев), столбик работает быстрее "Hьютона".
А также можно остановиться, получив нужное число знаков.
--
/* Alexandr Redchuk Re...@real.kiev.ua */
Andrey Vinokurov
Привет, Евгений.
В своем письме к "Valentin Shikhov" ты э..., как бы поделикатнее...,
выдал следующее:
VS> метод без деления на Х(i). Или побыстpее данного метода. А то деление
VS> такой медляк. :-(
EV>Пpи _точном_ вычислении коpня
"_Точно_" вычислить корень вообще невозможно, если аргумент не является
полным квадратом.
EV>избавиться от деления нельзя. Я думаю линейная
EV>аппpоксимация тебя не устpоит :).
Наивный чукотский юношь, ты нэ прав, аднака. Ты масква нэ биль, паравоз не
видэль, аднака.
Короче, не дезинформируй человека. _МОЖНО_ _ОБОЙТИСЬ_ _БЕЗ_ _ДЕЛЕНИЯ_.
Есть алгоритм извлечения кв. корня "столбиком" (похожий на деление
"столбиком") и его даже изучали в школе когда-то (мои родители, по крайней
мере, изучали его именно там). Этот алгоритм требует только
сравнений-вычитаний-сдвигов. И, что важно, его легко можно реализовать для
чисел повышенной точности. В своей простейшей (бинарной) форме он работает
с "родными", т.е. двоичными числами. Пару месяцев назад я закидывал сюда
описание этого алгоритма и мы дискутировали с рядом сомневающихся товарищей
по поводу того, "что лучше, Ньютон или столбик". Кстати, сама тема "isqrt"
- Integer SQuare RooT - осталась именно от тех времен.
Читай книжки, Женя. Они - кладезь разновсяческих знаний.
За сим усе.
Андрей.
Eugeny Vladimirtsev
Привет Andrey!
В четвеpг 09 Апpеля 11:27, Andrey Vinokurov писАл к All:
AV> Hаивный чукотский юношь, ты нэ прав, аднака. Ты масква нэ биль, паравоз не
AV> видэль, аднака.
Ж:-{
AV> Короче, не дезинформируй человека. _МОЖHО_ _ОБОЙТИСЬ_ _БЕЗ_
AV> _ДЕЛЕHИЯ_.
Знаю, знаю не кpичи так гpомко. Забыл.
AV> Читай книжки, Женя. Они - кладезь разновсяческих знаний.
И на стаpуху бывает...(с) Я сам году в 1982 ститал на "настольной pучной
вычислительной машине Феликс-1" квадpатные коpни этим методом Ж:-). У меня
где-то валяется инстpукция к аналогичной технике "ВК-1" с описанием этого
метода, но найти никак не могу.
С уважением, Eugeny
Alexandr A. Redchuck
Привет All и Винокуров Андрей!
Посмотрел я на $(SUBJ) столбиком. И как я его мог забыть :-)
(видать я учился по устаревшей программе, мне его таки в школе
давали, году так в 76-78).
Кто-то на днях писал, что это мол оптимизированный Hьютон, а
"оптимизировать можно все, что угодно". Hу не знаю. Мне так вылез
оптимизированный метод "Регистр Последовательных Приближений" aka SAR,
он же "метод поразрядного уравновешивания" (большинство микросхем
аналого-цифровых преобразователей на нем стоит). Из двоичного
варианта так и прет, в десятичном маскируется большим количеством
возможных весов разряда. А ньютона я в нем не нашел.
РПП:
поднимаем в регистре старший бит. Если функция содержимого регистра
меньше цели, оставляем этот бит, если больше - сбрасываем.
Повторяем вправо со всеми битами до младшего.
В тупом варианте для корня:
unsigned short sar_sqrt(unsigned long in) {
unsigned short sqr,mask=0x8000;
do {
sqr |= mask;
if( sqr*sqr > in ) sqr &= ~mask;
} while( mask >>= 1 );
return sqr;
}
Дальше для ухода от возведения в квадрат переходим к постепенному
приближению. Используем
(sqr+mask)*(sqr+mask) - sqr*sqr = 2*sqr*mask + mask*mask
и учитываем, что в mask поднята только одна единичка, т.е. умножение
на нее это сдвиг (в десятичном варианте таки надо множить -- стоИт
(10*2*R+t)*t и t приходится искать "на глаз" среди 10 весов разряда).
Можно гнать цикл по j от 15 до 0 и превратить
2*sqr*mask + mask*mask
в
(sqr << (j+1)) + ( 1 << (j+j) )
но при этом много сдвигов на переменную величину, которые не везде
легко делаются.
После оптимизации сдвигов
unsigned short isqrt( unsigned long in) {
unsigned long mask = 0x40000000, sqr = 0, temp;
do {
temp = sqr | mask;
sqr >>= 1;
if( temp <= in ) {
sqr |= mask;
in -= temp;
}
} while( mask >>= 2 );
// можно дать еще округление if( sqr > in ) ++sqr;
return (unsigned short)sqr;
}
"Школьный" вариант есть переходом к работе "цифра за цифрой", так как
при этом не нужно сравнивать младшие разряды. При этом вместо сдвига
маски разряда вправо двигается результат и вход влево.
"Цифра за цифрой" хороша тем, что можно получить любое число разрядов,
уйти вправо за точку. Hо для заранее известной и "небольшой" длины
может оказаться проще работать с полным числом, а не поразрядно (в
результате, а в источнике по паре разрядов), компьютеру-то все равно.
Вот мой вариант для i386
;------------------------------------------------------
; Integer square root calculation by ReAl
;------------------------------------------------------
; Parameter : unsigned long in stack
; Returns : unsigned short in AX
;------------------------------------------------------
.model tiny,C
TEXT SEGMENT PUBLIC 'CODE' PARA
ASSUME CS:TEXT
; .code
p386
isqrt proc C
public isqrt
mov edx,[esp+2] ; in
xor eax,eax ; sqr
bsr ecx,edx
jz done
and cl,1Eh
mov ebx,1
shl ebx,cl ; mask
loo:
lea ecx,[ebx+eax]
shr eax,1
cmp ecx,edx
jg short no
add eax,ebx
sub edx,ecx
no: shr ebx,2
jnz loo
done:
ret
isqrt endp
TEXT ENDS
end
По тактам выходит вроде бы быстрее, чем вариант "цифра за цифрой".
К Винокурову Андрею -- пожалуйста, проверь на своей тестовой
программе.
WBR,
Valery Datsjuk
Привет Victor!
В Среду 01 Апреля 1998, Victor Horbunkoff прислал письмо для Valentin Shikhov:
VH> Мнэ... Алгоритм деления в столбик чисел такого рода есть у дедушки
VH> Кнута. См. Д.Кнут "Искусство программирования для ЭВМ" т. 2
VH> "Получисленные алгоритмы" М., "Мир", 1977, стр.282-303. Корня я там,
VH> правда, не нашел. Кроме того: Ч.Уэзерелл "Этюды для программистов" М.,
VH> "Мир", 1982, стр.126-139 (умножение-деление, и ссылки на другие
VH> книги). Если тебя это интересовало, так вот оно.
А есть ли где Кнут в электронном виде?
С наилучшими пожеланиями, Well.
E-mail: we...@stack.ru IRC: SibNet #tomsk /Nick well