Názory k článku
Aritmetické operace s hodnotami ve formátu plovoucí řádové čárky

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)

Nebo Taylorovu radu pocitat pro pi/4 :-), tam sice vyjde dvakrat tolik nenulovych clenu, ale vypocet je na obe strany +- stejne presny.

Nechce se mi to přepočítávat, ale nevyjde náhodou Taylorova řada se středem v pi/2 víceméně stejně jako řada pro cosinus?

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.

"velmi nebezpečné a záludné A+B=A (platí samozřejmě pouze pro některé hodnoty A a B)"

To pro některé hodnoty A a B platí i v matematice. Teda spíš pro jednu konkrétní hodnotu B. (Samozřejmě mi je jasné, na co autor naráží.)

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).

0+\infty=NaN, zatimco \infty-\infty=\infty ?!

0+\infty je \infty a \infty-\infty je NaN, je to uvedene v prvni tabulce (druha kapitola).

Diky za odkaz. Sice se v tom clanku autor trosku navazi do me oblibene 6502, ale na druhou stranu vi o cem pise, takze na to ma pravo :-)

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

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)