Вычисление квадратного корня из целого числа
Николай Гарбуз
nick@sf.demos.su
Представленный алгоритм был создан в те бородатые времена, когда производительность x87 оставляла желать лучшего. Но и сейчас, скорость работы этого алгоритма соизмерима со скоростью вычисления с плавающей точкой на PII или MMX. На мой взгляд, материал может быть интересен, как начинающим программистам - пусть учатся писать программы, а не ломать их, и не вирусы, так и опытным - как игра ума., скомпилированный MSVC 5.0, на PII-233x2 дает следующие результаты (Листинг 1):
testing with range[0..1000] Done. testing range1=[0..1000]... fpu1... cpu1... cpu2... testing range2=[1000..10000]... fpu1... cpu1... cpu2... testing range3=[10000..100000]... fpu1... cpu1... cpu2... Done. range fpu1 cpu1 cpu2 1000 1.000 3.000 3.000 10000 1.000 3.000 4.000 100000 1.000 3.000 5.000
Задача вычисления квадратного корня при построении программ достаточна тривиальная. Функция для ее решения - sqrt - присутствует практически в любом из современных языков программирования. Однако практика использования функции sqrt показала, что данная функция ведет себя совершенно различным способом для целочисленных и действительных аргументов.
Пример 1
#include <stdio.h>
#include <math.h>
void main( )
{
int i = 169, j = 168;
printf(
<sqrt(%d)=%d, sqrt(%d)=%d>,
i, (int)sqrt(i),
j, (int)sqrt(j)
);
}
Результат выполнения кода приведенного в примере 1 выглядит так:
sqrt(169)=13, sqrt(168)=12
В действительности, значение квадратного корня для числа 168 соответствует числу 12.96, что по общепринятым правилам округления ближе к целому числу 13. В данном примере мы видим классический машинный случай округления с отбрасыванием дробной части.
Известно, многие сталкивались с подобной проблемой при построении целочисленных расчетных задач. Как правило, эта проблема решается написанием собственной функции, возвращающей правильный результат, т.е. результат, округленный до ближайшего целого, вместо отбрасывания дробной части. Один из вариантов собственной функции приведен в примере 2.
Пример 2
unsigned sqrt_fpu_true(long l)
{
unsigned rslt;
double f_rslt = 0;
if (l <= 0)
return 0;
f_rslt = sqrt(l);
rslt = (int)f_rslt;
if (!(f_rslt - rslt < .5))
rslt ++;
return rslt;
}
Функция, приведенная в примере 2, дает абсолютно правильные значения для всех целых чисел согласно принятым правилам округления. Однако возникает вопрос: возможно ли получение правильных результатов при использовании целочисленных алгоритмов?
Самый известный целочисленный алгоритм для вычисления квадратного корня из числа поражает своей простотой и приведен в примере 3.
Пример 3
unsigned sqrt_cpu_int(long l)
{ unsigned div = 1, rslt = 0;
while (l > 0)
{
l -= div, div += 2;
rslt += l < 0 ? 0 : 1;
}
return rslt;
}
Результат работы алгоритма из примера 3 идентичен результату из примера 1 - отбрасывание дробной части. Кроме того, невооруженным глазом виден еще один недостаток данного алгоритма - количество итераций в цикле соответствует значению вычисленного квадратного корня от аргумента L:
iteration count ~= sqrt(L) (1).
Рассматривая задачу вычисления квадратного корня с точки зрения уменьшения величины вычислительных затрат, более привлекателен целочисленный алгоритм, реализующий формулу Ньютона - пример 4.
Пример 4
unsigned sqrt_cpu_newton(long l)
{
unsigned rslt = (unsigned)l;
long div = l;
if (l <= 0)
return 0;
while (1)
{
div = (l / div + div) / 2;
if (rslt > div)
rslt = (unsigned)div;
else
return rslt;
}
}
Количество итераций в цикле для алгоритма из примера 4 приблизительно будет равняться натуральному логарифму от аргумента L:
iteration count ~= ln(L) (2).
Легко заключить, что разница в значениях формул (1) и (2) достаточно велика особенно для больших чисел, что и иллюстрирует ниже приведенная таблица.
Число L sqrt_cpu_int sqrt_cpu_newton 70000 264 11 300000 574 13 700000 836 13 990000 994 14
Однако результат работы алгоритма из примера 4 опять тот же - округление до целого числа отбрасыванием дробной части. Анализ кода алгоритма показывает, что наибольшая ошибка при вычислениях накапливается в главной формуле алгоритма и возникает при целочисленном делении на 2 без учета остатка от деления. В примере 5 приведен модифицированный алгоритм вычисления квадратного корня, с учетом вышеупомянутого замечания.
Пример 5
unsigned sqrt_cpu_newton(long l)
{ long temp, div = l;
unsigned rslt = (unsigned)l;
if (l <= 0)
return 0;
while (1)
{
temp = l / div + div;
div = temp >> 1;
div += temp & 1;
if (rslt > div)
rslt = (unsigned)div;
else
return rslt;
}
}
В модифицированный алгоритм добавлена одна переменная и две новые строки, реализующие целочисленное деление на 2 с учетом остатка. Модифицированный алгоритм вычисляет правильные значения - корень от аргумента с округлением до ближайшего целого практически для всех значений аргумента за исключением определенного ряда чисел. Для чисел этого ряда корень вычисляется, как число на единицу большее, чем истинное целочисленное его значение, определенное по общепринятым правилам округления (см. таблицу).
Число 2 6 12 20 30 42 56 72 90 110 132 Действит.корень 1,4 2,4 3,4 4,4 5,4 6,4 7,4 8,4 9,4 10,4 11,4 Целый корень 1 2 3 4 5 6 7 8 9 10 11 Вычисл.корень 2 3 4 5 6 7 8 9 10 11 12
Для всех аргументов алгоритма из приведенного ряда характерно, что вычисленное алгоритмом значение - есть целочисленный множитель, произведение которого с истинным целочисленным значением квадратного корня от аргумента дает сам аргумент. Знание выявленной закономерности позволяет сделать последнюю модификацию алгоритма, позволяющую окончательно устранить ошибки в вычислениях. Модификация заключается в проверке условия для выполнения коррекции результата вычисления при завершении работы алгоритма - пример 6.
Пример 6
unsigned sqrt_cpu_newton(long l)
{
long temp, div = l;
unsigned rslt = (unsigned)l;
if (l <= 0)
return 0;
while (1)
{
temp = l / div + div;
div = temp >> 1;
div += temp & 1;
if (rslt > div)
rslt = (unsigned)div;
else
{
if (l / rslt == rslt - 1 && l % rslt == 0)
rslt-;
return rslt;
}
}
}
Итак, вопрос о существовании целочисленного алгоритма для вычисления квадратного корня из целого числа с округлением результата до ближайшего целого по общепринятым правилам округления имеет утвердительный ответ.
В заключении следует отметить о существовании еще одной модификации алгоритма. На этот раз модификация преследует только задачу повышения производительности алгоритма. Повысить производительность итерационных алгоритмов возможно только одним способом - уменьшить количество итераций. Для приведенного в примере 6 алгоритма количество итераций можно значительно снизить, более точно подобрав начальные значения для переменной div - пример 7.
Пример 7
unsigned sqrt_newton(long l)
{
long temp , div;
unsigned rslt = (unsigned)l;
if (l <= 0)
return 0;
else if (l & 0xFFFF0000L)
if (l & 0xFF000000L)
div = 0x3FFF;
else
div = 0x3FF;
else
if (l & 0x0FF00L)
div = 0x3F;
else
div = (l > 4) ? 0x7 : l;
while (1)
{
temp = l / div + div;
div = temp >> 1;
div += temp & 1;
if (rslt > div)
rslt = (unsigned)div;
else
{
if (l / rslt == rslt - 1 && l % rslt == 0)
rslt-;
return rslt;
}
}
}
Последняя модификация алгоритма (пример 7) вычисляет квадратный корень из числа без ошибок округления на диапазоне [0..10000] в среднем за 3 итерационных цикла. В таблице ниже представлена сводная таблица по вычислительным затратам алгоритма на исследуемом диапазоне. На других диапазонах аргумента количество итераций не бывает больше 6, а в среднем равняется 3. Сравнивая с первоначально достигнутыми результатами, см. таблицу в начале, можно сказать, что достигнуто увеличение производительности как минимум в 2 - 4 раза.
Кол-во итераций 1 2 3 4 5 6 7 случаев из 10000 2 1965 6173 1779 80 0 0 % от всего 0,02% 19,65% 61,73% 17,19% 0,8% 0 0
Вероятно, что предел производительности алгоритма еще не достигнут, однако данная тема не является главной для настоящей статьи.
Листинг 1
// sqrt.cpp
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <conio.h>
unsigned sqrt_fpu1(long l)
{
if (l == 0)
return 0;
double f_rslt = sqrt(l);
unsigned rslt = (int)f_rslt;
if (!(f_rslt - rslt < .5))
rslt ++;
return rslt;
}
unsigned sqrt_cpu1(long l)
{
long temp;
unsigned div, rslt = l;
if (l <= 0)
return 0;
else if (l & 0xFFFF0000L)
if (l & 0xFF000000L)
div = 0x3FFF;
else
div = 0x3FF;
else
if (l & 0x0FF00L)
div = 0x3F;
else
div = (l > 4) ? 0x7 : l;
while (1)
{
temp = l / div + div;
div = temp >> 1;
div += temp & 1;
if (rslt > div)
rslt = div;
else
{
if (l / rslt == rslt - 1 && l % rslt == 0)
rslt--;
break;
}
}
return (unsigned)rslt;
}
unsigned sqrt_cpu2(long l)
{
if (l <= 0)
return 0;
long rslt = l, div = l;
while (1)
{
div = (l / div + div) / 2;
if (rslt > div)
rslt = div;
else
break;
}
return (unsigned)rslt;
}
unsigned sqrt_cpu3(long l)
{
unsigned div = 1;
unsigned rslt = 0;
while (l > 0)
{
l-= div, div += 2;
rslt += l < 0 ? 0 : 1;
}
return rslt;
}
unsigned sqrt_cpu4(long l)
{
unsigned div = 1, rslt = 0;
while (l > 0)
{
l-= div, div += 2;
rslt += l < 0 ? 0 : 1;
}
if (l != 0)
rslt++;
return rslt;
}
#define steps 1000
#define range1 1000L
#define range2 10000L
#define range3 100000L
#define count 2000
double times[3][3];
void CalcTime()
{
long l;
int i;
time_t first, second;
printf("testing range1=[%lu..%lu]...\n", 0L, range1);
printf("fpu1...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = 0l; l < range1; l += range1 / steps)
sqrt_fpu1(l);
second = time(NULL);
times[0][0] = difftime(second, first);
printf("cpu1...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = 0l; l < range1; l += range1 / steps)
sqrt_cpu1(l);
second = time(NULL);
times[0][1] = difftime(second, first);
printf("cpu2...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = 0l; l < range1; l += range1 / steps)
sqrt_cpu2(l);
second = time(NULL);
times[0][2] = difftime(second, first);
printf("testing range2=[%lu..%lu]...\n", range1, range2);
printf("fpu1...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = range1; l < range2; l += (range2 - range1) / steps)
sqrt_fpu1(l);
second = time(NULL);
times[1][0] = difftime(second, first);
printf("cpu1...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = range1; l < range2; l += (range2 - range1) / steps)
sqrt_cpu1(l);
second = time(NULL);
times[1][1] = difftime(second, first);
printf("cpu2...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = range1; l < range2; l += (range2 - range1) / steps)
sqrt_cpu2(l);
second = time(NULL);
times[1][2] = difftime(second, first);
printf("testing range3=[%lu..%lu]...\n", range2, range3);
printf("fpu1...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = range2; l < range3; l += (range3 - range2) / steps)
sqrt_fpu1(l);
second = time(NULL);
times[2][0] = difftime(second, first);
printf("cpu1...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = range2; l < range3; l += (range3 - range2) / steps)
sqrt_cpu1(l);
second = time(NULL);
times[2][1] = difftime(second, first);
printf("cpu2...\n");
first = time(NULL);
for (i = 0; i < count; i++)
for (l = range2; l < range3; l += (range3 - range2) / steps)
sqrt_cpu2(l);
second = time(NULL);
times[2][2] = difftime(second, first);
printf("Done.\n");
printf("range\t\t fpu1\t cpu1\t cpu2\n");
printf(
"%lu\t\t%5.3f\t%5.3f\t%5.3f\n",
range1,
times[0][0],
times[0][1],
times[0][2]
);
printf(
"%lu\t\t%5.3f\t%5.3f\t%5.3f\n",
range2,
times[1][0],
times[1][1],
times[1][2]
);
printf(
"%lu\t\t%5.3f\t%5.3f\t%5.3f\n",
range3,
times[2][0],
times[2][1],
times[2][2]
);
}
typedef unsigned (*sqrt_func)(long L);
void ViewDifferents(long rang1, long rang2, unsigned step, sqrt_func fpsqrt)
{
unsigned long l;
unsigned rf, ri;
printf("testing with range[%lu..%lu]\n", rang1, rang2);
for (l = rang1; l < rang2; l += (rang2 - rang1) / step)
{
rf = sqrt_fpu1(l);
ri = (*fpsqrt)(l);
if (rf != ri)
printf("sqrt(%lu) %u %u %6.2f\n", l, rf, ri, sqrt(l));
}
printf("Done.\n");
}
void main()
{
ViewDifferents(0, 1000, 1000, sqrt_cpu1);
CalcTime();
while (!kbhit());
}