Наши проекты:
Журнал · Discuz!ML · Wiki · DRKB · Помощь проекту |
||
ПРАВИЛА | FAQ | Помощь | Поиск | Участники | Календарь | Избранное | RSS |
[18.118.184.237] |
|
Сообщ.
#1
,
|
|
|
Может кто знает как организовать подобное. Нужно извлечь корень из числа, меняющегося от 0 до 1, причем единица это некая целая константа (0x7FFFFFFF например). Все должно быть в целых числах. Важна не столько точность алгоритма сколько его скорость. Скорость очень критична.
Всем заранее спасибо. |
Сообщ.
#2
,
|
|
|
народ!!!!!!!11
Подскажите хоть какой-нибудь алгоритм извлечения корня!!!!! |
Сообщ.
#3
,
|
|
|
Цитата 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 |
Сообщ.
#4
,
|
|
|
У меня есть две процедуры вычисления корня из 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 |
Сообщ.
#5
,
|
|
|
вот тебе несколько алг-мов
Цитата Вычисление квадратного корня из целого числа Николай Гарбуз 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> |
Сообщ.
#6
,
|
|
|
тебе тут стоко всего написали а самый нормальный алг. по системе ньютона с фиксированной точкой очень простой \\б берешь а(0)=1 в цикле а(т+1)= 0.5*(а(т)+1\а(т)) все т-номер иттерации сходимость очень быстрая квадратичная и не мудри с уважением Игорь Маленький |
Сообщ.
#7
,
|
|
|
2 Z@ если критична скорость, то без таблицы начальных приближений при любом итерационном процессе не обойтись. Максимальная скорость достигается тогда, когда в памяти хранится вся таблица значений. (Известная дилемма - память супротив скорости ). Построй таблицу начальных приближений. Само число - индекс в таблице. Быстро находишь первое приближение с нужной точностью, а дальше - формула Ньютона. В своё время на ЕС ЭВМ скорость вычисления sqrt повысилась на ~40\%
Кроме того, неплохо бы знать под какой процессор нужна оптимизация? Задержка / пропускная способность (в тактах) для FSQRT (DP) на Intel Pentium 4 =38/38. Так что ускорять нечего. |