 |
Математика: Быстрое вычисление функций и констант.
Вычисление квадратного корня из целого числа
Автор: Николай Гарбуз,
nick@sf.demos.su
Представленный алгоритм был создан в те бородатые времена, когда производительность x87 оставляла желать лучшего. Но и сейчас, скорость работы этого алгоритма соизмерима со скоростью вычисления с плавающей точкой на PII или MMX. На мой взгляд, материал может быть интересен, как начинающим программистам - пусть учатся писать программы, а не ломать их, и не вирусы, так и опытным - как игра ума. Тест, скомпилированный MSVC 5.0, на PII-233x2 дает следующие результаты:
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
#include
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 (l)
{
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 (l)
{
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 (l)
{
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) reslt--;
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 (l)
{
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) reslt--;
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 |
Вероятно, что предел производительности алгоритма еще не достигнут, однако
данная тема не является главной для настоящей статьи.
Вверх по странице, к оглавлению и навигации.