5

Подскажите, пожалуйста, алгоритм нахождения квадратного корня из числа (не sqrt(a) и не pow(a, 0.5)), а итерациями, вручную. Ответ должен быть не целочисленным, а точным, типа float.

Язык - желательно C.

1
  • в гугле нет алгоритмов, дающих целочисленный ответ. а здесь, возможно, кто-то сталкивался. Но за гостеприимство и отзывчивость в любом случае спасибо =/
    – Kollibry
    8 апр 2012 в 15:33

3 ответа 3

20

Можно попробовать метод Ньютона. По сути, чтобы найти корень из числа S, надо решить уравнение f(x) = 0, где

Подставляя эту функцию в формулу для итераций метода Ньютона, получаем

Я набросал работающую программу на Python, можете по ней ориентироваться:

#!/usr/bin/python
import sys

s = float(sys.argv[1])

x = 1.0
while abs(x * x - s) > 0.00001:
    x = (x * x + s) / 2. / x  
print x
3
  • Спасибо! буду разбираться
    – Kollibry
    8 апр 2012 в 15:51
  • еще можно с помощью ряда тэйлора(см. вики) 9 апр 2012 в 16:31
  • 2
    Классическая запись: x=0.5(x+S/x). 21 ноя 2015 в 1:28
9

Есть бывший школьный алгоритм извлечения квадратного корня столбиком, по типу деления уголком. Сохранён для потомков в Википедии, п.6.6. Алгоритм построен на тождестве
(ae+b)2=a2e2 + (2ae+b)b для e=10.

В роли a выступает значение квадратного корня, вычисленное для старших цифр числа, в роли b - новая цифра. Значение корня вычисляется по недостатку, количество десятичных цифр результата не ограничено.

В отличие от алгоритма деления столбиком, в получении новой цифры участвуют сразу две новых цифры исходного числа, а при получении дробной части нули добавляют тоже парами. Об этом нам говорит множитель e2=100 в формуле.

Чтобы результат не получил множитель sqrt(10)=3.162, цифры попарно группируют, начиная от десятичной точки. При этом в самой старшей группе оказываются одна или две цифры.

введите сюда описание изображения На рисунке (пример из Википедии, рисунок изменён) изображён процесс извлечения квадратного корня из числа N = 69696.

Число N размещено привычным образом. Вместо уголка справа мы писали знак =. И почти сразу писали первую цифру b2=a2=2, потому что [sqrt(6)]=2. Это был разряд сотен.

На второй итерации рассуждаем так: левая часть тождества - сверху, число a2=2 - справа от знака равенства. Вычитая из верхнего числа величину (a2)2=4 и снося к этой двойке вниз пару цифр 96, мы получаем соотношение (2a2e+b1)b1 = 296, где b1 - неизвестная цифра в разряде десятков. Проводим горизонтальную и вертикальную чёрточки, как на рисунке, и слева от вертикальной черточки с отступом влево на один разряд пишем число 2a2=4.
Если дописать справа цифру b1, получим как раз 2a2e+b1. А если эту же цифру b1 подписать внизу, то получим пример на умножение с ответом 296. Подберём эту цифру в уме: 45х5=225 - маловато, 47х7=329 - уже перебор, 46х6=276 - в самый раз! И сразу записываем 276 под 296 - чтобы не забыть. А праведно добытую цифру b1=6 дописываем к ответу справа.

Итак, третья итерация. И кое-что уже ясно: две новые чёрточки, вычитание столбиком и снос следующих двух цифр - 2096... А вот и нюанс: (2a2e+b1)+b1=2(a2+b1)=2a1. Т.е. мы не обязаны умножать в уме 26х2, а можем сложить столбиком 46 и 6. И получить 52.

Итак, 52b0 x b0 = 2096... b0=4 - ура, сошлось! Аккуратно подрисовываем 4 в умножение, 2096 снизу, 4 справа. Корень извлечён в целости и сохранности нацело.

Теперь ничто не мешает формализовать алгоритм, причём для любой системы счисления.

Может возникнуть вопрос: зачем этот алгоритм, если есть итерационный алгоритм Герона Александрийского? Просто потому, что у него есть плюсы - такие же, как у деления уголком:

  1. Возможность вычисления результата со скоростью ~10 операций (сравнение, вычитание, логическое ИЛИ и сдвиги) на бит результата с аппаратной погрешностью.
  2. Удобство применения для длинных двоичных чисел, представленных массивами.
  3. Возможность извлечения квадратного корня вручную в любой системе счисления.
  4. Гарантированное представление результата по недостатку.

Например, для e=2 программа эмулирует аппаратную реализацию квадратного корня для мантиссы числа (программа написана на PHP):

function msb($num){
    $k = ($num > 65535) ? 16 : 0;   // Если левое слово ненулевое, это минимум 16 разрядов.
    if (($num>>$k) > 255) $k |= 8;  // Второй байт после сдвига не 0 - это +8 разрядов.
    if (($num>>$k) > 15)  $k |= 4;  // Левая тетрада после сдвига не 0 - это +4 разряда.
    if (($num>>$k) > 3)   $k |= 2;  // Левая диада после сдвига не 0 - это +2 разряда
    if (($num>>$k) > 1)   $k |= 1;  // Левый разряд после сдвига не 0 - это +1 разряд.
    return $k;
}

function sqrt_school($n){
    if($n>>31) return null; 
    if($n<1) return 0;
    $digit_pos = msb($n) >> 1;              // текущая позиция бита результата 
    $result = 1 << $digit_pos;              // результат
//  printf ("n= %b digit_pos=%d", $n, $digit_pos); ////
    $top = $n - ($result<<$digit_pos);      // остаток
    while($digit_pos>0){
        $mult = ($result << $digit_pos--) | (1 << ($digit_pos<<1)); // вычитаемое для цифры 1
//      printf ("<br>result=%b digit_pos=%d top=%b mult=%b", $result, $digit_pos, $top, $mult); 
        if($top >= $mult){
            $result |= 1 << $digit_pos;     // вписали бит в результат
            $top -= $mult;                  // уменьшили остаток
        }
//      printf("&emsp;left=%b",$left);
    }
    return $result; 
}

function sqrt_heron($n){
    $s = (float)$n;
    $x = (float)1;
    while (abs($x*$x - $s) > 0.5) {
        $x = ($x + $s/$x)/2;
    }  
    return (int)$x;
}

for($i=0; $i<10000; $i++){
    $test_array[$i] = mt_rand(2,1000);  
} 
$time0 = microtime(GET_AS_FLOAT);
for($i=0; $i<10000; $i++){
    $school = sqrt_school($test_array[$i]);
} 
$time1 = microtime(GET_AS_FLOAT);
for($i=0; $i<10000; $i++){
    $heron = sqrt_heron($test_array[$i]);
} 
$time2 = microtime(GET_AS_FLOAT);
printf("<br>Test [2,1000]: &emsp; time_school=%d ns &emsp; time_heron=%d ns",100*($time1-$time0)+0.5,100*($time2-$time1)+0.5);

for($i=0; $i<10000; $i++){
    $test_array[$i] = 1000000*mt_rand(2,1000);  
} 
$time0 = microtime(GET_AS_FLOAT);
for($i=0; $i<10000; $i++){
    $school = sqrt_school($test_array[$i]);
} 
$time1 = microtime(GET_AS_FLOAT);
for($i=0; $i<10000; $i++){
    $heron = sqrt_heron($test_array[$i]);
} 
$time2 = microtime(GET_AS_FLOAT);
printf("<br>Test 1000000000*[2,1000]: &emsp; time_school=%d ns &emsp; time_heron=%d ns",100*($time1-$time0)+0.5,100*($time2-$time1)+0.5);

for($i=1; $i<3000000000; $i*=3){
    $school = sqrt_school($i);
    $heron = sqrt_heron($i);
    printf("<br>sqrt_school($i)=%d &emsp; sqrt_heron($i)=%d" , $school, $heron);    
} 

Попутно проведено сравнение с алгоритмом Герона (замена деления на двойку умножением ничего не даёт).
Результаты совпали, по быстродействию в 3-5 раз выигрывает "школьный" алгоритм.

Test [2,1000]:   time_school=13 ns   time_heron=39 ns
Test 1000000 * [2,1000]:   time_school=19 ns   time_heron=90 ns
sqrt_school(1)=1   sqrt_heron(1)=1
sqrt_school(3)=1   sqrt_heron(3)=1
sqrt_school(9)=3   sqrt_heron(9)=3
sqrt_school(27)=5   sqrt_heron(27)=5
sqrt_school(81)=9   sqrt_heron(81)=9
sqrt_school(243)=15   sqrt_heron(243)=15
sqrt_school(729)=27   sqrt_heron(729)=27
sqrt_school(2187)=46   sqrt_heron(2187)=46
sqrt_school(6561)=81   sqrt_heron(6561)=81
sqrt_school(19683)=140   sqrt_heron(19683)=140
sqrt_school(59049)=243   sqrt_heron(59049)=243
sqrt_school(177147)=420   sqrt_heron(177147)=420
sqrt_school(531441)=729   sqrt_heron(531441)=729
sqrt_school(1594323)=1262   sqrt_heron(1594323)=1262
sqrt_school(4782969)=2187   sqrt_heron(4782969)=2187
sqrt_school(14348907)=3787   sqrt_heron(14348907)=3787
sqrt_school(43046721)=6561   sqrt_heron(43046721)=6561
sqrt_school(129140163)=11363   sqrt_heron(129140163)=11363
sqrt_school(387420489)=19683   sqrt_heron(387420489)=19683
sqrt_school(1162261467)=34091   sqrt_heron(1162261467)=34091

Из тестов следует также, что увеличение разрядности результата по "школьному" алгоритму не приведёт к серьзному замедлению вычислений.

Частичный перевод описания алгоритма на английский язык здесь

2
  • Стоит привести код, либо пример, т.к. ссылка даже на википедию может внезапно оказаться недоступной.
    – insolor
    25 ноя 2015 в 16:00
  • 1
    @insolor Получите - распишитесь! :) 25 ноя 2015 в 19:18
1

Да, Вам нужен не целочисленный. Если кому-то интересен целочисленный, то вот (вроде хорощий) из интернета

unsigned int 
sqrt32(unsigned long n)  
{  
    unsigned int c = 0x8000;  
    unsigned int g = 0x8000;

    for(;;) {  
        if(g*g > n)  
            g ^= c;  
        c >>= 1;  
        if(c == 0)  
            return g;  
        g |= c;  
    }  
}
1
  • Я не уверен, но IMHO если этот код модифицировать для 64-bit, а потом применить к мантиссе из double (скорректировав ее и порядок для нечетного порядка (т.к. порядок делим на 2)), а после вычисления нормализовать и опять упаковать мантиссу и порядок в double, то должна вычислять. Но, проверять надо. А вот смысл задания я не особо понимаю. Видимо учебное для схемотехников или моделирование аппаратуры ? В железе делается что-то в этом духе.
    – avp
    9 апр 2012 в 8:42

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge you have read our privacy policy.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.