Подскажите, пожалуйста, алгоритм нахождения квадратного корня из числа (не sqrt(a)
и не pow(a, 0.5)
), а итерациями, вручную. Ответ должен быть не целочисленным, а точным, типа float
.
Язык - желательно C.
Подскажите, пожалуйста, алгоритм нахождения квадратного корня из числа (не sqrt(a)
и не pow(a, 0.5)
), а итерациями, вручную. Ответ должен быть не целочисленным, а точным, типа float
.
Язык - желательно C.
Можно попробовать метод Ньютона. По сути, чтобы найти корень из числа 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
Есть бывший школьный алгоритм извлечения квадратного корня столбиком, по типу деления уголком. Сохранён для потомков в Википедии, п.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 справа. Корень извлечён в целости и сохранности нацело.
Теперь ничто не мешает формализовать алгоритм, причём для любой системы счисления.
Может возникнуть вопрос: зачем этот алгоритм, если есть итерационный алгоритм Герона Александрийского? Просто потому, что у него есть плюсы - такие же, как у деления уголком:
Например, для 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(" 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]:   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("<br>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("<br>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
Из тестов следует также, что увеличение разрядности результата по "школьному" алгоритму не приведёт к серьзному замедлению вычислений.
Частичный перевод описания алгоритма на английский язык здесь
Да, Вам нужен не целочисленный. Если кому-то интересен целочисленный, то вот (вроде хорощий) из интернета
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;
}
}