Nebylo, protoze cos() se pocita taky Taylorovou radou. Autor ukazal nevhodny postup - natvrdo 3 cleny rozvoje. V praxi je v compute_sin cyklus, ktery probiha tak dlouho, dokud velikost dalsiho clenu rozvoje presahne pozadovanou presnost. Navic je dobre si uvedomit, ze nasledujici clen rozvoje se snadno da spocist z predchoziho. Vubec neni nutne pocitat n-tou mocninu a faktorial n.
Ten demonstracni priklad na vypocet sinu nebyl myslen jako ukazka presne implementace v FPU, pouze to je demonstrace, ze i s pouhymi tremi cleny Taylorovy rady se da dopidit k rozumne presnosti. V FPU je, pokud vim (ted nemluvim o x86), vetsinou pocet clenu take omezen, to zavisi od pouzivaneho formatu (single, double).
Faktorialy se nepocitaji vubec, pro tech par clenu jsou primo ulozeny v FPU (stejne, jako hodnoty uhlu pro CORDIC) a mocniny se samozrejme daji prevest na opakovane nasobeni.
Samozřejmě, že bylo, protože pro x z leveho okoli pi/2 je potřeba na požadovanou přesnost více členů Taylorova rozvoje. Takže je rozumné numericky počítat Taylorovy řady jen pro x z intervalu (0,pi/4)
V pi/2 ano, ja jsem psal pi/4, coz je polovina intervalu, na kterem se sin() pocita. To (pokud to trosku zjednodusim) znamena, ze chyby budou od pi/4 na obe strany rust, ale maximalni chyba (ta je u sin() i cos() na konci intervalu) bude cca polovicni.
Prednosti pouziti pi/4 je dale to, ze se nemusi pocitat tak velke mocniny a faktorialy (ale ty jsou stejne v tabulce), protoze zadny ze clenu neni nulovy, narozdil od sin(a=0) a cos(a=0). Nevyhodou jsou ty Pisvajcovy konstanty diky tomu, ze sin(pi/4), cos(pi/4) atd. nejsou moc "zaokrouhlena" cisla.
Vlastne jo, nejak jsem zapomel, ze i pro realy plati A+B=A, ovsem jen pro mnozinu nulove miry (u FP je tech B-cek, pro ktere vztah plati, bohuzel prilis mnoho).
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);
}
To 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.
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)