sqrt с фиксированной точкой

Может кто знает как организовать подобное. Нужно извлечь корень из числа, меняющегося от 0 до 1, причем единица это некая целая константа (0x7FFFFFFF например). Все должно быть в целых числах. Важна не столько точность алгоритма сколько его скорость. Скорость очень критична.
Всем заранее спасибо.

народ!!!!!!!11
Подскажите хоть какой-нибудь алгоритм извлечения корня!!!!!

Цитата Z@, 30.07.02, 01:26:12

Подскажите хоть какой-нибудь алгоритм извлечения корня!!!!!

type <br> TNum = 0..1;<br> var  a: TNum <br> implementation<br>procedyre Sqrt1; <br>var<br>b:real;<br>Begin<br>a:=StrToFloat(Edit1.Text);<br>b:=Sqrt(a);<br>Edit1.Text:=FloatToStr(b);<br>end;<br>

Любой так любой ;D

У меня есть две процедуры вычисления корня из fixedpoint 16:16 (не проверял).


;--------------------------------------------------
; SqrtLP - Fixed Point Square Root (Low Precision)
;    IN     : ecx
;   OUT     : eax
;  Modified : ebx,ecx,edx

SqrtLP           PROC
;NOTE:
;this function only gives a 8.8 precision !!!!!!!!!!!!!

     xor     eax,eax

     mov     ebx,40000000h
sqrtLP1:
     mov     edx,ecx         ;edx = val
     sub     edx,ebx         ;val - bitsqr
     jl      sqrtLP2
     sub     edx,eax         ;val - root
     jl      sqrtLP2
     mov     ecx,edx         ;val >= (root+bitsqr) -> accept subs
     shr     eax,1           ;root >> 1
     or      eax,ebx         ;root | bitsqr
     shr     ebx,2           ;bitsqr>>2
     jnz     sqrtLP1
     shl     eax,8
     ret
sqrtLP2:
     shr     eax,1           ;val < (root+bitsqr) -> dont change val
     shr     ebx,2           ;bitsqr>>2
     jnz     sqrtLP1
     shl     eax,8
     ret

SqrtLP           ENDP




;--------------------------------------------------
; Sqrt - Fixed Point Square Root (High/Normal Precision)
;    IN     : ecx
;   OUT     : eax
;  Modified : ebx,ecx,edx
Sqrt         PROC

;This is the High Precision version for the sqrt.
;It gives the optimal 8.16 precision but takes
;a little longer (24 iterations, 48 bits intead of
;16 iterations and 32 bits)

     xor     eax,eax         ;eax is root
     mov     ebx,40000000h
sqrt1:
     mov     edx,ecx         ;edx = val
     sub     edx,ebx         ;val - bitsqr
     jb      sqrt2
     sub     edx,eax         ;val - root
     jb      sqrt2
     mov     ecx,edx         ;val >= (root+bitsqr) -> accept subs
     shr     eax,1           ;root >> 1
     or      eax,ebx         ;root | bitsqr
     shr     ebx,2           ;bitsqr>>2
     jnz     sqrt1
     jz      sqrt5
sqrt2:
     shr     eax,1           ;val < (root+bitsqr) -> dont change val
     shr     ebx,2           ;bitsqr>>2
     jnz     sqrt1
; we now have the 8.8 precision

sqrt5:
     mov     ebx,00004000h
     shl     eax,16
     shl     ecx,16
sqrt3:
     mov     edx,ecx         ;edx = val
     sub     edx,ebx         ;val - bitsqr
     jb      sqrt4
     sub     edx,eax         ;val - root
     jb      sqrt4
     mov     ecx,edx         ;val >= (root+bitsqr) -> accept subs
     shr     eax,1           ;root >> 1
     or      eax,ebx         ;root | bitsqr
     shr     ebx,2           ;bitsqr>>2
     jnz     sqrt3
     ret
sqrt4:
     shr     eax,1           ;val < (root+bitsqr) -> dont change val
     shr     ebx,2           ;bitsqr>>2
     jnz     sqrt3
     ret

Sqrt           ENDP

вот тебе несколько алг-мов

Цитата

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

Николай Гарбуз
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());
}

Цитата

Вариант вычисления квадратного корня (алгоритм Ньютона)
Для вычисления квадратного корня здесь использован известный метод Ньютона. Рассмотрены машинно-зависимый и машинно-независимый варианты. Перед применением алгоритма Ньютона, область определения исходного числа сужается до [0.5,2] (в другом варианте до [1,16]). Второй вариант машинно-независим, но работает дольше.

Скажем сразу, не самый быстрый вариант, но один из самых быстрых. Основная проблема заключается здесь в том, что нет универсального машинного представления чисел с плавающей точкой, поэтому разделение числа на мантиссу и двоичную экспоненту как составляющих его компьютерного представления нельзя записать единым образом. Имено поэтому здесь подключено описание математической библиотеки, из которой, впрочем, используются только frexp() и ldexp(). Конкретная реализация этих функций очень проста, но машинно-зависима. При разделении числа на мантиссу и экспоненту, мантисса оказывается в пределах [0.5,1). Если при этом экспонента нечетная, мантисса умножается на 2, тем самым область определения становится [0.5,2). К мантиссе применяется алгоритм Ньютона с начальным приближением 1, и затем окончательный результат получается с помощью применения ldexp() с половинным значением экспоненты.

Для машинно-независимого варианта выделение экспоненты заменяется серией последовательных делений исходного значения на 16, пока аргумент не станет определяться на интервале [1,16]. Если исходный аргумент был меньше 1, то перед серией делений обращаем его. Алгоритм Ньютона применяется с начальным приближением 2. Окончательный результат получается серией последовательных умножений на 4 и дальнейшим обращением, если аргумент обращался.

Сам алгоритм Ньютона для вычисления a=Sqroot(x) представляет быстро сходящуюся (при хорошем начальном приближении) серию итераций: ai+1=0.5*(ai+x/ai), где i -- номер итерации.

Ниже приводятся программы, реализующие вычисление квадратного корня по указанным алгоритмам.


--------------------------------------------------------------------------------

/* Square root almost without mathematic library. Optimized for floating
  point single precision. The functions frexp() and ldexp() from the
  mathematic library are very fast and machine-dependent.
  The second variant is less fast but does not use the math library
  at all.
  Copyright © Nikitin V.F. 2000

  Square root decomposing the argument into the mantisse and the
  exponent (faster program):
  float Sqroot(float x);

  Square root without usage of any external functions (slower but
  machine-independent):
  float Sqroot1(float x);

  In case of invalid domain (x<0.) functions return zero (do not
  generate an error as library equivalents).
*/

#include <math.h>   /* frexp() and ldexp() only */

/* 4 iterations needed for the single precision */
#define ITNUM 4

/* variant using external f.p. decomposition/combining to [0.5,1] */
float Sqroot(float x) {
 int expo,i;
 float a,b;
 if(x<=0.F) return(0.F);
 /* decompose x into mantisse ranged [0.5,1) and exponent. Machine-
    dependent operation is presented here as a function call. */
 x=frexp(x,&expo);
 /* odd exponent: multiply mantisse by 2 and decrease exponent
    making it even. Now the mantisse is ranged [0.5,2.) */
 if(expo&1) {x*=2.F;expo--;}
 /* initial approximation */
 a=1.F;
 /* process ITNUM Newtonian iterations with the mantisse as
    an argument */
 for(i=ITNUM;i>0;i--) {
   b=x/a; a+=b; a*=0.5F;
 }
 /* divide the exponent by 2 and combine the result.
    Function ldexp() is opposite to frexp. */
 a=ldexp(a,expo/2);
 return(a);
}

/* variant without math lib usage. The domain is shrunk to [1,16].
  Repeated divisions by 16 are used. */
float Sqroot1(float x) {
 int sp=0,i,inv=0;
 float a,b;
 if(x<=0.F) return(0.F);
 /* argument less than 1 : invert it */
 if(x<1.F) {x=1.F/x;inv=1;}
 /* process series of division by 16 until argument is <16 */
 while(x>16.F) {sp++;x/=16.F;}
 /* initial approximation */
 a=2.F;
 /* Newtonian algorithm */
 for(i=ITNUM;i>0;i--) {
   b=x/a; a+=b; a*=0.5F;
 }
 /* multiply result by 4 : as much times as divisions by 16 took place */
 while(sp>0) {sp--;a*=4.F;}
 /* invert result for inverted argument */
 if(inv) a=1.F/a;
 return(a);
}

мой (какой-то из этих)

<br>CLNum sqrt(const CLNum& x)<br>{<br>      int counter=0,i;<br>      bool inverted=false;<br><br>      CLNum nolik("0");<br>      CLNum c16((char)16);<br>      CLNum c05("0.5");<br>      CLNum c4((char)4);<br>      CLNum a(x);<br>      if(a<nolik) <br>            throw(CLNum::Error("MathError: Value that is less than zero is under sqr"));<br>      if(a==nolik)return a;//корень из нуля.<br>      /* argument less than 1 : invert it */<br>      if(a<(CLNum((char)1))) <br>      {<br>            a=((CLNum)(char)1)/a;<br>            inverted=true;<br>      }<br>      <br>      /* process series of division by 16 until argument is <16 */<br>      while(a>c16) <br>      {<br>            counter++;<br>            a/=c16;<br>      }<br>      /* initial approximation */<br>      CLNum p1((char)2);<br>      CLNum p2(1,0);<br>      /* Newtonian algorithm */<br>      for(i=ITNUM;i>0;i--) <br>      {<br>//            cout<<"p2="<<p2<<"  p1="<<p1<<"  a="<<a<<"\n";<br>//          cin>>r11;<br>            p2=a/p1; <br>//            cout<<"p2="<<p2<<"\n";<br>//            cin>>r11;<br>            p1+=p2; <br>//            cout<<"p1="<<p1<<"\n";<br>//            cin>>r11;<br>            p1*=c05;<br>      }<br>      /* multiply result by 4 : as much times as divisions by 16 took place */<br>      while(counter>0) <br>      {<br>            counter--;<br>            p1*=c4;<br>      }<br>      /* invert result for inverted argument */<br>      if(inverted) p1=((CLNum)(char)1)/p1;<br>      return(p1);<br><br>}<br>

тебе тут стоко всего написали а самый нормальный алг. по системе ньютона с фиксированной точкой очень простой \\б

берешь а(0)=1
в цикле а(т+1)= 0.5*(а(т)+1\а(т))


все т-номер иттерации

сходимость очень быстрая квадратичная




и не мудри

с уважением Игорь Маленький

2 Z@ если критична скорость, то без таблицы начальных приближений при любом итерационном процессе  не обойтись. Максимальная скорость достигается тогда, когда в памяти хранится вся таблица значений. (Известная дилемма - память супротив скорости ). Построй таблицу начальных приближений. Само число - индекс в таблице. Быстро находишь первое приближение с нужной точностью, а дальше - формула Ньютона. В своё время на ЕС ЭВМ скорость вычисления sqrt повысилась на ~40\%
Кроме того, неплохо бы знать под какой процессор нужна оптимизация?
Задержка / пропускная способность (в тактах) для FSQRT (DP) на Intel Pentium 4 =38/38.
Так что ускорять нечего.