Vlákno názorů k článku Aritmetické operace s hodnotami ve formátu plovoucí řádové čárky od BruXy - Jelikoz 1.0/2.0*(xi+y/xi) ma tu blbou vlastnost, ze se...
-
BruXy (neregistrovaný)Jelikoz 1.0/2.0*(xi+y/xi) ma tu blbou vlastnost, ze se tam musi delit, a deleni jako nebrat, da se vyuzit pro vypocet odmocniny reseni rovnice 1/y*y - A = 0 (Newtonovou metodou), potom je iteracni vztah y_{i+1} = 1.0/2.0*y_i*(3-A*y_i*y_i), pomoci ktereho se da urcit hodnota 1/sqrt(A) a pak sqrt(A) = A * 1/sqrt(A).
Presnost vypoctu zavisi na volbe y_0, ta musi byt v rozmezi 0 < y_0 < sqrt(3/a) a vlastne celkovy pocet iteraci k dosazeni presneho vysledku zavisi na volbe pocatecniho odhadu (muzeme ho trefit od oka vypoctem nebo mit odhady predpocitane v tabulce).
#include <stdio.h>
#include <math.h>
double sqrt_step(float y, float A){
return 1.0/2.0*y*(3-A*y*y);
}
int main(){
double vstup = 10000;
double yi = ((3/vstup) / 2); // odhad ma byt 0 < y0 < sqrt(3/vstup);
int i;
double sqr=1.0/sqrt(vstup); // počitadlo iterac
for( i = 0; i < 20; i++){
yi = sqrt_step(yi, vstup);
printf("%d\t%11.5f\t%11.5f\t%10.2f%%\n",
i+1,
yi, // vypočtená hodnota
fabs(yi - sqr), // absolutní chyba
100.0f * fabs(yi-sqr) /sqr ); // relativní chyba }
}
printf("sqrt(%11.5f) = %11.5f\n", vstup, vstup * yi);
}
1 0.00022 0.00978 97.75%
2 0.00034 0.00966 96.63%
3 0.00051 0.00949 94.94%
4 0.00076 0.00924 92.42%
5 0.00114 0.00886 88.65%
6 0.00170 0.00830 83.05%
7 0.00252 0.00748 74.81%
8 0.00370 0.00630 63.02%
9 0.00529 0.00471 47.05%
10 0.00720 0.00280 28.00%
11 0.00893 0.00107 10.66%
12 0.00984 0.00016 1.65%
13 0.01000 0.00000 0.04%
14 0.01000 0.00000 0.00%
15 0.01000 0.00000 0.00%
16 0.01000 0.00000 0.00%
17 0.01000 0.00000 0.00%
18 0.01000 0.00000 0.00%
19 0.01000 0.00000 0.00%
20 0.01000 0.00000 0.00%
sqrt(10000.00000) = 100.00000 -
Pavel TišnovskýZlatý podporovatelTo vypada zajimave, diky za uverejneni. Jen se chci zeptat - neexistuje tam vice reseni te rovnice? (jeste jsem tu funkci 1/y*y-A nezkoumal, ale mam dojem, ze ma aspon dva koreny). Chyba - mam aspon ten dojem - klesa trosku pomaleji, ale ta eliminace deleni IMHO stoji za to. -
BruXy (neregistrovaný)Rekl bych, ze zavisi na volbe pocatecniho y_0. Ten rozsah 0 < y_0 < sqrt(3/A) je vybran obecne tak, ze
to vyjde kladne doufam pro vsechna A > 0. Ale vyzkousej dosadit jiny pocatecni odhad, treba 0.005 pro A = 10000, spravny vysledek po 5 iteracich. O duvod vic pouzit predpocitanou tabulku (klidne muze byt i v ROMce daneho FPU)
