Очень нужна помощь в оптимизации алгоритма извлечения квадратного корня Опции

[Rostislav]

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

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

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 циклов.

Сообщение отредактировал Rostislav - Dec 5 2012, 07:35

[Guest_TSerg_*]

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

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

Еще раз приведу (для 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;

[Rostislav]

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

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

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

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

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

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

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

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

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

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

Сообщение отредактировал Rostislav - Dec 21 2012, 11:16

[beaRTS]

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

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

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

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

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

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

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

Код

    while( firstEstimation - nextEstimation != 0 ) {    
        firstEstimation = nextEstimation;

        nextEstimation = ( firstEstimation + quotient )/2;

        //for next iteration check in while()
        quotient = radicand/firstEstimation;

        //for debug
        ++iterationCnt;
    } //while

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

Код

    for( int k = 0, g = 50; k < 50; ++k, --g) {
        temp = k*k + g*g;
        i = ( float )k;
        q = ( float )g;
        std::cout << "<cmath> sqrt(): " << sqrt(temp) << "     my sqrt(): " << sqrtIQ(i, q) << "   perform in " << iterationCnt << " iterations"<< std::endl;
    } // for

Эскизы прикрепленных изображений