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

[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

[alexeyv]

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

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

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

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

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

[Rostislav]

1 = 1^2
1+3 = 2^2
1+3+5 = 3^2
и так далее.

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

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

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

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

Так?

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

[Tanya]

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

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

[Rostislav]

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

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

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

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

[MrYuran]

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

Код

........b=b+(b+2)

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

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

[Rostislav]

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

Код

........b=b+(b+2)

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

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

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

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

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

[skyv]

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

[Rostislav]

Я использовал просто таблицу

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

[Lmx2315]

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

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

[MrYuran]

 bibliofond_556017.rtf ( 117.89 килобайт ) Кол-во скачиваний: 1040


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

[Rostislav]

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

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

[sysel]

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

Код

// Вычисление квадратного корня из 32х битного числа
int16 isqrt32 (int32 from){
  register int32 mask=0x40000000;
  register int32 sqr=0;
  do {
    register int32 temp = sqr | mask;
    sqr >>= 1;
    if (temp <= from){
      sqr |= mask;
      from -= temp;
    };
  } while (mask >>= 2);
  if (sqr < from) sqr++;
  return ((int16)sqr);
}

[skyv]

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

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

[анатолий]

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

Код

function SQROOT(x1:integer) return integer is
        variable m,n,i, a, x,y,b: INTEGER;      
        constant bn:integer:=28;
            begin                
        y:=x1;     x:=x1;
        m:=1;    n:=bn/2;
        L1: for i in 1 to n loop
            a:=4*x;     
            exit L1 when a>=2**(bn);    
            y:=2*y;    x:=a;
        end loop;    
        L2:for i in 1 to n-1 loop
            b:=x+x/2**m;
            a:=b+b/2**m;
            if a<2**(bn) then
                x:=a;
                y:=y+y/2**m;
            end if;    
            m:=m+1;
        end loop;        
        return y/(2**n);
    end;

[Rostislav]

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

Код

....
    b:=x+x/2**m;
....

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

[Guest_TSerg_*]

Код

function intSqrt_(Val : Longword) : Longword;
var bitSqr : Longword;
begin
bitSqr := $40000000;
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;

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

Код

function Q(x: Longword): Longword;
var k: Longword;
begin
  Result := 0;
  k := 1;
  while (k <= x) do begin
    Dec(x,k);
    Inc(k,2);
    Inc(Result);
  end;
end;

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

[akorud]

Есть возможность использовать 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интезируется в логическую схему.

[VitekSVM]

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

[ASN]

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

[Rostislav]

С какой практической точностью нужно вычислять выражение ?

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

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

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

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

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

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

----

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

Код

  a : Word;  // подкоренное выражение
  b : Word;  // результат
  n : Byte;

  a := 121;
  b := 0;

  // начало
  for n := 15 downto 0 do   // Счет от 15 до 0   // 15 - разрядность подкоренного выражения (вернее N не 15, а N=15+1)
   begin
     if ((b + 1 shl n) * (b + 1 shl n) < a) or         // shl - побитовый сдвиг влево на n бит
        ((b + 1 shl n) * (b + 1 shl n) = a) then
      begin
        b := b + 1 shl n;
      end;
   end;
  // конец

  Form1.Caption := IntToStr (b);    // b = 11

  Циклов для 16-ти разрядного числа 16!!!
         для 32-х  разрядного числа 32!!!!!!

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

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

[beaRTS]

[code] a : Word; // подкоренное выражение
b : Word; // результат
n : Byte;

a := 121;
b := 0;

// начало
for n := 15 downto 0 do // Счет от 15 до 0 // 15 - разрядность подкоренного выражения (вернее N не 15, а N=15+1)
begin
if ((b + 1 shl n) * (b + 1 shl n) < a) or // shl - побитовый сдвиг влево на n бит
((b + 1 shl n) * (b + 1 shl n) = a) then
begin
b := b + 1 shl n;
end;
end;
// конец

Form1.Caption := IntToStr (; // b = 11
[code]

По-моему, получилось красиво)))

Может, я чего-то не понимаю... но в цикле 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;
// конец

Сообщение отредактировал beaRTS - Dec 21 2012, 02:10

[Rostislav]

temp := b + 1 shl n;
if (temp * temp < a) or //А ТАК НЕЛЬЗЯ????

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

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

[beaRTS]

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

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

[Rostislav]

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

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

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

[blackfin]

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

[Rostislav]

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

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

[beaRTS]

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

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

[Rostislav]

...я эту ссылку первый раз вижу...

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

[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

[Rostislav]

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

?

Сообщение отредактировал Rostislav - Dec 22 2012, 08:27

[beaRTS]

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

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

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

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

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

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

[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

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

 

[beaRTS]

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

[Guest_TSerg_*]

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

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

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

bitSqr := $FFFF;

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