Квадратный корень из числа

#php #cpp #алгоритм #математика #числа

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

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

Ответы

Ответ 1

Можно попробовать метод Ньютона. По сути, чтобы найти корень из числа 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


Ответ 2

Есть бывший школьный алгоритм извлечения квадратного корня столбиком, по типу деления
уголком. Сохранён для потомков в Википедии, п.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 справа. Корень извлечён в целости и сохранности нацело.

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

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


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


Например, для 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 ("
result=%b digit_pos=%d top=%b mult=%b", $result, $digit_pos,
$top, $mult); 
        if($top >= $mult){
            $result |= 1 << $digit_pos;     // вписали бит в результат
            $top -= $mult;                  // уменьшили остаток
        }
//      printf(" 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("
Test [2,1000]:   time_school=%d ns   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("
Test 1000000000*[2,1000]:   time_school=%d ns   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("
sqrt_school($i)=%d   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


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


Ответ 3

Да, Вам нужен не целочисленный. Если кому-то интересен целочисленный, то вот (вроде
хорощий) из интернета
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;  
    }  
}