поиск целочисленного квадратного корня

Тема в разделе "WASM.A&O", создана пользователем Mikl___, 25 окт 2010.

  1. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.708
    Казалось бы, что может быть проще fsqrt, но мы пойдем другим путем
    1) Самый известный целочисленный алгоритм для вычисления квадратного корня из числа поражает своей простотой
    unsigned sqrt_cpu_int(usigned X)
    { unsigned div = 1, result = 0;
    while (X > 0)
    { X -= div; div += 2; result += X < 0 ? 0 : 1; }
    return result;
    }
    Недостаток данного алгоритма - количество итераций будет увеличиваться с ростом Х
    2) вычисление квадратного корня по разложению в ряд Тейлора.
    Пусть Х - любое число; f(X) - некоторая функция, зависящая от X; a - известное число, близкое к Х; f(a) - известное значение функции.
    Разложим f(X) в ряд Тейлора:
    f(X)=f(a)+(X-a)f '(a)+((X-a)² *f" (a))/2! + ... + ((X-a)n *f n (a))/n!
    Пусть X - число, из которого нужно извлечь корень. Тогда f(X)=√[o]X[/o]; a - известное число близкое к X; f(a)=√[o]a[/o] - известное число близкое к √[o]X[/o], и f(X)=√[o]a[/o] +(X-a)/(2√[o]a[/o])+...=(2a+X-a)/(2√[o]a[/o])=(a+X)/(2√[o]a[/o])
    Величина √[o]X[/o] может быть найдена, если задаться величиной √[o]a[/o] и затем вычислить f(X). f(X)² можно сравнить с исходным числом Х. Если точность окажется недостаточной, тогда число а заменяется на f(X)², √[o]a[/o] на f(X) и вычисление повторяется
    3) поиск целочисленного квадратного корня методом Ньютона начинается с некоторого значения g0, которое является начальной оценкой √[o]X[/o]. Затем выполняется серия уточнений значения квадратного корня по формуле gn+1=(gn+X/gn)/2. Для уменьшения количество итераций можно на первом этапе более точно подобрать начальное значения для переменной g0
    { usigned x1;
    int a, g0, g1;
    if ( x <= 1 ) return x;
    a = 1;
    x1 = x - 1;
    if (x1 > 65535) {a = a + 8; x1 = x1 >> 16;}
    if (x1 > 255) {a = a + 4; x1 = x1 >> 8;}
    if (x1 > 15) {a = a + 2; x1 = x1 >> 4;}
    if (x1 > 3) {a = a + 1;}
    g0 = 1 << a; // g0 = 2ª
    g1 = (g0 + (x >> a)) >> 1;
    while (g1 < g0) // повторяем, пока приближения строго уменьшаются
    { g0 = g1; g1 = (g0 + (x/g0)) >> 1}
    }
    4) Так как деление достаточно медленная операция, можно от нее отказаться. Для вычисления квадратного корня используем умножение. 32-разрядное число Х будет иметь 16-разрядный квадратный корень. На каждом шаге происходит уточнение 1 бита значения корня. Для ускорения сделано "безбранчевое" вычисление прибавить или отнять значение в соответствующем разряде
    Код (Text):
    1.     mov ebx,4000h
    2.     mov ecx,8000h
    3.     mov eax,40000000h
    4. @@: cmp eax,X
    5.     je @f; преждевременный выход из цикла
    6.     sbb edx,edx; если CF=1 тогда edx=0FFFFFFFFh иначе edx=0
    7.     sub ecx,ebx
    8.     and edx,ebx; если CF=1 тогда edx=ebx иначе edx=0
    9.     lea ecx,[ecx+edx*2]; если CF=1 тогда ecx=ecx+ebx иначе ecx=ecx-ebx
    10.     shr ebx,1; переходим к следующему разряду
    11.     mov eax,ecx
    12.     mul eax; получаем "eax в квадрате"
    13.     test ebx,ebx
    14.     jnz @b
    15. @@: mov eax,ecx; eax:=sqrt(X)
     
  2. Y_Mur

    Y_Mur Active Member

    Публикаций:
    0
    Регистрация:
    6 сен 2006
    Сообщения:
    2.494
    Mikl___
    Меня тоже дико восхитил способ 1 :)) Думаю есть способы его улучшения.
    По своей сути этот метод - способ пробежаться по ряду квадратов целых чисел 1, 4, 9, 16, 25, 36... и имхо делать эту пробежку можно более эффективно чем вычитанием.
     
  3. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.708
    Y_Mur
    [​IMG]
    Это сумма арифметической прогрессии из нечетных чисел, а как вам 4-ый способ?
    [offtop]рисунок сгенерировался на сайте http://dxdy.ru/math/ по строкам class="latex" alt="$\sum\limits_{n=1}^{k}(2n-1)=k^2$" title="$\sum\limits_{n=1}^{k}(2n-1)=k^2$" /> Вот бы и на WASM.RU/FORUM такое же![/offtop]
     
  4. Y_Mur

    Y_Mur Active Member

    Публикаций:
    0
    Регистрация:
    6 сен 2006
    Сообщения:
    2.494
    Mikl___
    Не сразу вник в 4) - да метод отличный :)
     
  5. gorodon

    gorodon New Member

    Публикаций:
    0
    Регистрация:
    19 окт 2009
    Сообщения:
    301
    Mikl___
    Метод 4 хороший... суть - частный метод решения уравнения y=f(x) на отрезке [x2,x1] методом деления отрезка пополам.
    Пример кода (рекурсивное решение):
    Код (Text):
    1. //x1 > x2
    2. long sqrt_long(long x_2, long x1, long x2)
    3. {
    4.     if(x1 == x2 + 1)
    5.         return x1;
    6.     //
    7.     long x = (x1+x2)>>1;//(x1+x2)/2
    8.     //
    9.     if(x*x >= x_2) return sqrt_long(x_2, x, x2);
    10.     //
    11.     return sqrt_long(x_2, x1, x);
    12. }
    13. ...
    14.     long x_2 = ...
    15.     long x = sqrt_long(x_2,65536,0);
    16.     printf("sqrt(%d)=%d",x_2,x);
     
  6. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    Я так понимаю, 3) и 4) (c) Уоррен? Там и ещё один метод есть, без умножений/делений, как и 1), но логарифмической сложности - функции isqrt4/isqrt5 в http://www.hackersdelight.org/HDcode/isqrt.c.txt .
     
  7. qqwe

    qqwe New Member

    Публикаций:
    0
    Регистрация:
    2 янв 2009
    Сообщения:
    2.914
    gorodon
    много итераций.

    можно попробовать ограничить диапазон учитывая, что старший бит X в X^2 = Y, где значащие биты Y в диапазоне [N, 0], будет на N/2 позиции. те подбирать надо биты только на участке [N/2 - 1, 0]

    АДД
    а лучше даже не половинным делением, а тупым побитным перебором с проверкой двигаясь от N/2 - 1 к 0

    на 64 бита [63, 0] это будет максимум 30 итераций

    (а может и нет. не проверял)

    ЗЫ извинямс, если ктото уже


    АДД2:
    Код (Text):
    1. unsigned sqrt(unsigned Y){
    2.   unsigned N = get_last_bit_pos(Y);
    3.   unsigned X = 1 << (N >> 1);
    4.   unsigned bit;
    5.  
    6.   if(N == 0) return Y;
    7.  
    8.   for(bit = X >> 1; bit > 0; bit >>= 1){
    9.     X |= bit;
    10.     if(X * X > Y) X &= ~bit;
    11. // тут можно еще на равенство проверять
    12. // if(X * X == Y) return X;
    13.   }
    14.  
    15.   return X;
    16. }
     
  8. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.708
    3) Да, Уоррен ("Алгоритмические трюки для программистов" -- М. :Издательский дом "Вильямс" 2003 стр 198 Листинг 11.1 Метод Ньютона поиска целочисленного квадратного корня) 4) мой, но не скрою, написал после прочтения Уоррена "Аппаратная реализация поиска целочисленного квадратного корня" и "Бинарный поиск целочисленного квадратного корня"
     
  9. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.708
    Если уж речь зашла об Уоррене, то вот метод который меня "убил"
    то же на ассемблере
    Код (Text):
    1.     mov ebp,X
    2.     mov ecx,40000000h
    3.     xor eax,eax
    4. @@: mov ebx,eax
    5.     or ebx,ecx
    6.     shr eax,1
    7.     cmp ebp,ebx
    8.     sbb edx,edx
    9.     not edx
    10.     mov edi,edx
    11.     and edx,ebx
    12.     and edi,ecx
    13.     sub ebp,edx
    14.     or eax,edi
    15.     shr ecx,2
    16.     test ecx,ecx
    17.     jnz @b
    18. ;eax=sqrt(X)
    Не понятно как, ни деления, ни умножения, но работает! Причем алгоритм был описан аж в 1945 в Von Neumann J. "First Draft of a Reaport on the EDVAC"
     
  10. gorodon

    gorodon New Member

    Публикаций:
    0
    Регистрация:
    19 окт 2009
    Сообщения:
    301
    Дык это все тот же алгоритм 4 - побитовой обработки... :)
    Квадрат от каждого бита записан в m (10^2=100 100^2=10000... (bin)) - потом нужные биты "съезжают" на нужное число разрядов вниз, уже находясь в y.
    Это как дифференцирование реализовывать с помощью интеграторов... :) - из той же серии оптимизации...
     
  11. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.708
    Чуть исправил код #9
    Код (Text):
    1. mov ebp,X
    2.     inc ebp
    3.     mov ecx,40000000h; здесь для уменьшения количества итераций
    4. ;можно подбирать значение, как в методе Ньютона
    5.     xor eax,eax
    6. @@: lea ebx,[eax+ecx]
    7.     shr eax,1
    8.     cmp ebx,ebp
    9.     sbb edx,edx
    10.     mov edi,edx
    11.     and edx,ebx
    12.     sub ebp,edx
    13.     and edi,ecx
    14.     or eax,edi
    15.     shr ecx,2
    16.     jnz @b
     
  12. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.708
    Код (Text):
    1.     mov ebp,X
    2.         inc ebp
    3.     bsr ecx,ebp
    4.     and ecx,-2
    5.     mov ebx,1
    6.     shl ebx,cl;для уменьшения количества итераций
    7.         xor eax,eax
    8. @@: lea ecx,[eax+ebx]
    9.     shr eax,1
    10.     cmp ecx,ebp
    11.     sbb edx,edx
    12.     mov edi,edx
    13.     and edx,ecx
    14.     sub ebp,edx
    15.     and edi,ebx
    16.     or eax,edi
    17.     shr ebx,2
    18.     jnz @b