Jump to content

    
Sign in to follow this  
Rostislav78

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

Recommended Posts

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;

Share this post


Link to post
Share on other sites
Топик-стартер слушает только себя в попытке изобрести велосипед.:)

 

Вам уже приводились идеи и работающие варианты.

Науке известны разные способы извлечения целочисленного корня на разные вкусы.

 

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

 

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

 

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

 

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

 

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

 

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

 

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

 

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

Edited by Rostislav

Share this post


Link to post
Share on other sites

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

 

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

 

          1  2. 3  4 
        ___________
       /
     \/  01 52.27 56

         01                   1*1 <= 1 < 2*2                 x = 1
         01                     y = x*x = 1*1 = 1
         00 52                22*2 <= 52 < 23*3              x = 2
         00 44                  y = (20+x)*x = 22*2 = 44
            08 27             243*3 <= 827 < 244*4           x = 3
            07 29               y = (240+x)*x = 243*3 = 729
               98 56          2464*4 <= 9856 < 2465*5        x = 4
               98 56            y = (2460+x)*x = 2464*4 = 9856
               00 00          Algorithm terminates: Answer is 12.34

          1. 4  1  4  2
          _____________
       /
     \/  02.00 00 00 00

         02                  1*1 <= 2 < 2*2                 x = 1
         01                    y = x*x = 1*1 = 1
         01 00               24*4 <= 100 < 25*5             x = 4
         00 96                 y = (20+x)*x = 24*4 = 96
            04 00            281*1 <= 400 < 282*2           x = 1
            02 81              y = (280+x)*x = 281*1 = 281
            01 19 00         2824*4 <= 11900 < 2825*5       x = 4
            01 12 96           y = (2820+x)*x = 2824*4 = 11296
               06 04 00      28282*2 <= 60400 < 28283*3     x = 2
                             The desired precision is achieved:
                             The square root of 2 is about 1.4142

?

Edited by Rostislav

Share this post


Link to post
Share on other sites
ЗЫ. Кстати, очень интересно какие все же идеи есть у beaRTS)))

да какие там идеи... =) тож начитался много чего (по диагонали)

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

 

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

 

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

 

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

 

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

Share this post


Link to post
Share on other sites
ЗЫ. Кстати, очень интересно какие все же идеи есть у 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

post-69111-1356166744_thumb.png

Share this post


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

 

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

 

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

 

bitSqr := $FFFF;

 

Да, чисто механически оттранслировал на 16р:)

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

 

 

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this