Перейти на главную   
  helloworld.ru - документация и книги по программированию  
helloworld.ru - документация и книги по программированию    
    главная     хостинг     создание сайтов    
Поиск по сайту:  
Смотрите также
Языки программирования
C#
MS Visual C++
Borland C++
C++ Builder
Visual Basic
Quick Basic
Turbo Pascal
Delphi
JavaScript
Java
PHP
Perl
Assembler
AutoLisp
Fortran
Python
1C

Интернет-технологии
HTML
VRML
HTTP
CGI
FTP
Proxy
DNS
протоколы TCP/IP
Apache

Web-дизайн
HTML
Дизайн
VRML
PhotoShop
Cookie
CGI
SSI
CSS
ASP
PHP
Perl

Программирование игр
DirectDraw
DirectSound
Direct3D
OpenGL
3D-графика
Графика под DOS

Алгоритмы
Численные методы
Обработка данных

Системное программирование
Драйверы

Базы данных
MySQL
SQL

Другое

Хостинг


Друзья
demaker.ru
Реклама

Лучший хостинг. Аренда серверов




helloworld.ru

Вычисление квадратного корня из целого числа

Николай Гарбуз

Представленный алгоритм был создан в те бородатые времена, когда производительность 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());
}










helloworld.ru © 2001-2016
Все права защищены
Rambler's Top100 TopList Rambler's Top100