Ребята, всем привет! Помогите, пожалуйста оптимизировать алгоритм извлечения корня. А то что-то до фига получается циклов вычисления. Заранее огромное спасибо всем откликнувшимся!

Алгоритм простой, основан на методе последовательного приближения. Собственно вот он (без шелухи):

bprev:=bprev div 2;
if (b+bprev)*(b+bprev) > a then
b:=b-bprev
else
if (b+bprev)*(b+bprev) < a then
b:=b+bprev
else
выход;
повтор цикла....

Полный текст:

a:word;
b:word;
bprev:word;
cnt:word;

a:=strtoint(edit1.Text); // подкоренное выражение
b:=a; // в b хранится результат
bprev:=a;
cnt:=0;
while 1=1 do begin
cnt:=cnt+1; // Счетчик циклов
bprev:=bprev div 2;
if bprev=0 then begin
bprev:=1;
if ((b+1)*(b+1) > a) and ((b+0)*(b+0) < a) then break;
end;
if (b+bprev)*(b+bprev) > a then
b:=b-bprev
else
if (b+bprev)*(b+bprev) < a then
b:=b+bprev
else
break;
end;
b:=b+bprev;
Form1.Caption:='результат: '+inttostr(b)+' циклов: '+inttostr(cnt);

Результаты замера следующие:
sqr (1) - 2 цикла;
sqr (10) - 5 циклов;
sqr (100) - 5 циклов;
sqr (1000) - 13 циклов;
sqr (10000) - 38 циклов;
sqr (32000) - 73 цикла;
...

Алгоритм предложенный alexeyv:

i:=1; // результат в i
b:=1;
a:=10000; // подкоренное выражение
while 1=1 do begin
if (b > a) or (b = a) then break;
b:=b+2;
a:=a-b;
i:=i+1;
end;
caption:=inttostr(i);

sqr (10000) - 100 циклов.

Если нужна только целая часть квадратного корня, то можно воспользоваться другим алгоритмом:

Для квадратов чисел верны следующие равенства:
1 = 1^2
1+3 = 2^2
1+3+5 = 3^2
и так далее.

То есть, узнать целую часть квадратного корня числа можно, вычитая из него все нечётные числа по порядку, пока остаток не станет меньше следующего вычитаемого числа или равен нулю, и посчитав количество выполненных действий. Например, так:
Выполнено 3 действия, квадратный корень числа 9 равен 3.

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

Циклов конечно будет больше, зато время их выполнения гораздо меньше

Забавный алгоритм))

Т.е. если я Вас правильно понял, то программа должна выглядеть так:

a=подкоренное выражение
b=1
i=1

цикл
если b > a или b = a то
....выход (результат в i)
иначе
....b=b+(b+2)
....i=i+1
повтор цикла

Так?

Что Вас так забавляет?
Немного школьной арифметики...
Сумма арифметической прогрессии дает - сумма N нечетных чисел в точности равна N в квадрате.
Следующий разряд (после точки) можно вычислить как (половина остатка)/N.

Что забавляет? Все)) Магия чисел в особенности))

Проверил работу алгоритма:

sqr (10000) = 13
Но sqr (9) = 3 !!!!!

Где-то ошибка(

i:=1;
b:=1;
a:=10000;
while 1=1 do begin
if (b > a) or (b = a) then break;
b:=b+(b+2);
i:=i+1;
end;
caption:=inttostr(i);

просто b=b+2
a=a-b

(эх и ужос этот ваш паскаль.. да ещё без тега code)

Кстати, детская арифметика заканчивается начиная с sqr (36). Результат 5!



Ага, спасибо, все заработало! В итоге sqr (10000) вычисляется за 100 шагов. Оптимизации не получилось. Можно попробовать изменить его под переменную базу. Сейчас мы прибавляем 2, а можно начинать, например, с 16 последовательно уменьшая это число в двое. Если конечно это сработает. Буду думать...



Ну да, вот так и мучаюсь всю жизнь!))

Я использовал просто таблицу, которую предварительно вычислял.
Где ее лучше расположить это зависит от ресурсов CPU.

А Вы можете привести пример? Мне нужно как-то оценить эффективность такого подхода. Сейчас я не могу сделать никаких выводов. Я знаю, что существуют табличные вычисления, но все нужно проверять в живую)

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

Нажмите для просмотра прикрепленного файла

Можно скомбинировать.
Таблицу для грубого определения и дальше уточнять.

Спасибо! Изучаю Вавилонский способ. Оказывается древние не только воевать умели) Я пока с комбинированием не хочу связываться. Нужно выжать максимум из того, что сейчас наработано. Мне кажется, что потенциал есть у обоих алгоритмов. А о комбинировании различных методов мысля уже проскальзывала, но хочется сделать монолитно-красиво)

Не знаю, как называется метод, но пашет:

Так и есть.
В этом случае операция вычисления (в принципе чего угодно) сводится к
чтению ячейки массива.

Вот функция корня по алгоритму CORDIC

Анатолий, спасибо огромное! Изучаю алогоритмик. Только вот я не понял, что означает в программе '**'? Степень числа?

Если сравнивать с уже упомянутой медленной функцией:


то разница на писюках более чем в 20 раз лучше.

Есть возможность использовать floating-point? Может известный алгоритм быстрого инверсного корня из Quake III http://en.wikipedia.org/wiki/Fast_inverse_square_root и факт что sqrt(x) = x/sqrt(x) -> sqrt(x) = x * rsqrt(x)

Это степень числа, язык VHDL. Cинтезируется в логическую схему.

http://algolist.manual.ru/maths/count_fast/index.php

Rostislav
С какой практической точностью нужно вычислять выражение ?
На какой целевой платформе выполняется прикладная программа?
Если вычислительная система поддерживает вычисления с плавающей точкой , то лучше воспользоваться арифметическим сопроцессором (типа сопроцессора для i386) или специальными командами самого процессора (как в TMS320C55).
Если требуется целочисленная арифметика, то проще и быстрее заранее вычисленной таблицы в памяти что-то предложить сложно.

Уважаемый ASN! Мне придумалась еще одна идея. Более простого алгоритма я придумать не смог (кстати, не понятно, почему этот вариант мне раньше не пришел в голову). Но все же отвечу на Ваши вопросы))

Меня интересует целочисленный результат.



Меня интересовал вариант, который мог бы выполнятся на микроконтроллерах с простой архитектурой (типа PIC12).



Таблицы сразу не понравились. Для них нужна память.

----

Итак, мое заднее)))) предложение:



По-моему, получилось красиво))) Операции умножения можно заменить на бинарное сложение столбиком с предварительным смещением бит влево. Если более оптимального алгоритма наука)) еще не придумала, то тему можно закрывать. Теперь возьмусь за вычисление тригонометрических функций)))) Всем спасибо!!!!!

Может, я чего-то не понимаю... но в цикле for явно не красота красуется. Вы 5 раз вычисляете одно и тоже выражение b + 1 shl n !! Дык, зачем это делать?. один раз вычислите, присвойте результат какой-нибудь временной переменной и таскайте ее за собой везде. что-то типа:

for n := 15 downto 0 do // Счет от 15 до 0 // 15 - разрядность подкоренного выражения (вернее N не 15, а N=15+1)
begin
temp := b + 1 shl n;
if (temp * temp < a) or // А ТАК НЕЛЬЗЯ???? (temp * temp <= a )
(temp * temp = a) then
begin
b := temp
end;
end;
// конец

Да нее, Вы все правильно понимаете)))))))) Просто я не стал делать очевидное, то, что лежит на поверхности. Ну конечно же, мне известны некоторые трюки программирования. Умножение на 2 (и деление), например, легко можно заменить на побитовый сдвиг и т.п. А Ваш вариант оптимизации сделает, скорее всего, компилятор. Меня (и не только меня, наверное) больше всего интересовала суть. Суть я и изложил. Поэтому, спасибо Вам за участие и с наступающим Новым Годом!!!!!!!)))))))))

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

Здорово!!!))) Выкладывайте Ваши идеи!!!!!! Очень интересно!!!)))

Сто раз уже обсуждали: ReAl.

Да, похоже, я не открыл Америку))))))))

повторение - мать учения.. к тому же поиск хренова-то работает на форуме. я эту ссылку первый раз вижу. Так что спасибо за "пароли, явки, адреса").

Присоединяюсь! К тому же, может быть "подует свежий ветер"!)))

Топик-стартер слушает только себя в попытке изобрести велосипед.

Вам уже приводились идеи и работающие варианты.
Науке известны разные способы ивлечения целочисленного корня на разные вкусы.

Еще раз приведу (для 16р сетки), число иттераций всегда 9 и никаких умножений.

function QuickSqrt_(Val : word) : word;
var bitSqr : word;
begin
bitSqr := $10000;
Result := 0;

While bitSqr <> 0 do begin
if (Val >= bitSqr + Result) then begin
Dec(Val, bitSqr + Result);
Result := (Result shr 1) or bitSqr;
end
else Result := Result shr 1;
bitSqr := bitSqr shr 2;
end;
end;

Уважаемый TSerg, ну тё Вы ругаетесь? Себя я, конечно же, слушаю))) Велосипед? Да, изобретаю! Не мешайте мне побыть гениальным

Ваш вариант действительно очень эффективен. Да, выполняется за всего за 9 циклов! Я бы добавил его в начало темы, но не могу. Первое сообщение не доступно уже для изменения.

Более того, в последнем Вашем коде есть ошибка: bitSqr : word. Это говорит о том, что Вы лично его не проверяли в деле.

Судя по тому, что маске присваивается значение 0x010000, маска bitSqr должна быть 32-х разрядной (как и было в первом Вашем варианте). Далее в коде все операции с ней происходят как с 32-х разрядным числом (компилятор не будет динамически менять разрядность), это требует дополнительной памяти данных и инструкций. Это первое.

Второе. Объем кода в Вашем варианте мне показался больше, чем в моем (к тому же, учитывая операции с 32-х разрядным числом в микроконтроллере с 8-ми разрядной архитектурой (даже с аппаратным умножителем 8х8), например). Но это на вскидку. Если все совсем делать серьезно, то, по идее, оба алгоритма нужно перевести на ассемблер микроконтроллера PIC16, например, и сравнить.

Но, Ваш код прост. Мой имеет один (?) недостаток - операция умножения. В микроконтроллере с аппаратным умножителем это не будет проблемой, в остальных да.

Вариантов, действительно, было множество. Но ищется оптимальный (имеющий малый объем кода, простоту вычислений, быстрый). Не все они отвечают этим требованиям. И вкус, часто определяется архитектурой вычислительного ядра. Ищу золотую середину)))

ЗЫ. Кстати, очень интересно какие все же идеи есть у beaRTS)))

TSerg, Ваш вариант описывается тут: http://en.wikipedia.org/wiki/Methods_of_co...ng_square_roots

Ребят, кто-нибудь понимает предложенные там еще два метода:

?

да какие там идеи... =) тож начитался много чего (по диагонали)
мне корень нужен для нахождения модуля вектора, т.е.: имею синфазную и квадратурную составляющую сигнала (I и Q), т.е. как бы координаты вектора, ну и надо модуль найти, чтобы два канала I и Q соединить в один. Так вот. У меня на проце TMS320F28335 есть FPU (в общем, floating-point понимает и перемножитель даже есть аппаратный), работаю с числами float (вещественными одинарной точности).

Ну и показался метод Ньютона (он же вроде вырождается в метод Вавилонский) более подходящий из всех, т.к. ему пофигу что считать: либо целые либо вещественны числа. Стал думать как делать для него первое приближение g0 (в книге Уоррен Генри. "Алгоритмические трюки для программистов такие обозначения" ). Ну а тут прям в масть ознакомился с неким алгоритмом "максимум альфа плюс минимум бета", который линейно аппроксимирует выражение sqrt( I^2 + Q^2) с приемлемой точностью и быстротой (точность зависит от коэффициентов альфа и бета). ну и решил первое приближение g0 делать при помощи него.

итого, программка работает. только мне еще вкурить надо как должным образом проверить ее производительность во времени (что-то типа профилирования), сколько итераций расчета по методу ньютона производит с таким начальным/первым приближением. и сравнить по времени выполнения со стандартными функциями С++. Ну и величину ошибки тоже вычислить.

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

В общем статистику еще не проводил и не оформлял. но стопудово сделаю какую-нить статейку.

да какие там идеи... =) тож начитался много чего (по диагонали)
мне корень нужен для нахождения модуля вектора, т.е.: имею синфазную и квадратурную составляющую сигнала (I и Q), т.е. как бы координаты вектора, ну и надо модуль найти, чтобы два канала I и Q соединить в один. Так вот. У меня на проце TMS320F28335 есть FPU (в общем, floating-point понимает и перемножитель даже есть аппаратный), работаю с числами float (вещественными одинарной точности).

Ну и показался метод Ньютона (он же вроде вырождается в метод Вавилонский) более подходящий из всех, т.к. ему пофигу что считать: либо целые либо вещественны числа. Стал думать как делать для него первое приближение g0 (в книге Уоррен Генри. "Алгоритмические трюки для программистов такие обозначения" ). Ну а тут прям в масть ознакомился с неким алгоритмом "максимум альфа плюс минимум бета", который линейно аппроксимирует выражение sqrt( I^2 + Q^2) с приемлемой точностью и быстротой (точность зависит от коэффициентов альфа и бета). ну и решил первое приближение g0 делать при помощи него.

итого, программка работает. только мне еще вкурить надо как должным образом проверить ее производительность во времени (что-то типа профилирования), сколько итераций расчета по методу ньютона производит с таким начальным/первым приближением. и сравнить по времени выполнения со стандартными функциями С++. Ну и величину ошибки тоже вычислить.

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

В общем статистику еще не проводил и не оформлял. но стопудово сделаю какую-нить статейку.

ну вот сырые данные в картинке во вложении.
Под итерациями я здесь подразумеваю сколько раз алгоритму требуется пройтись по циклу while() :



а в main() прогоняю такие числа:

по поводу своей идеи начал новую тему. вот здесь кому интересно - зайдите , посоветуйте

Это не мой вариант, но давно известный.

bitSqr := $FFFF;

Да, чисто механически оттранслировал на 16р
Delphi исправляет мою погрешность, но в микроконтроллере надо быть внимательнее, конечно же.