C/Zaawansowane operacje matematyczne: Różnice pomiędzy wersjami
Usunięta treść Dodana treść
Nie podano opisu zmian |
m →Obliczenia numeryczne: kolejność |
||
Linia 118:
{{Uwaga|Należy brać pod uwagę, że w komputerze liczby rzeczywiste nie są tym samym, czym w matematyce. Komputery nie potrafią przechować każdej liczby zmiennoprzecinkowej, w związku z tym obliczenia prowadzone przy użyciu komputera mogą być niedokładne i odbiegać od prawidłowych wyników.<ref>[http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html|What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg]</ref> Szczególnie ważne jest to przy programowaniu aplikacji inżynieryjnych oraz w medycynie, gdzie takie błędy mogą skutkować katastrofą i/lub narażeniem ludzkiego życia oraz zdrowia.<ref>[http://www.ima.umn.edu/~arnold/455.f97/notes.html Two disasters caused by computer arithmetic errors]</ref>}}
===Epsilon maszynowy===
Linia 954 ⟶ 586:
lnew = 19000000000 new = 1820130816
</pre>
=====Zapobieganie=====
Linia 1175 ⟶ 812:
double frac = r - (int)r ;
</syntaxhighlight>
===Błędy w obliczeniach numerycznych===
Typy:<ref>[https://www.sci.brooklyn.cuny.edu/~mate/nml/numanal.pdf INTRODUCTION TO NUMERICAL ANALYSIS WITH C PROGRAMS by Attila Mate]</ref><ref>[https://people.eecs.berkeley.edu/~demmel/cs267/lecture21/lecture21.html Basic Issues in Floating Point Arithmetic and Error Analysis by Jim Demmel ]</ref>
* wg etapu operacji :<ref>[https://www.securecoding.cert.org/confluence/display/c/FLP32-C.+Prevent+or+detect+domain+and+range+errors+in+math+functions securecoding.cert.org : Prevent+or+detect+domain+and+range+errors+in+math+functions]</ref>
** blędne dane wejściowe : niezgodne z oczekiwanym typem
** dane wejściowe powodują błąd rezultatu
* wg rodzaju operacji
** zmiennoprzecinkowych<ref>[http://stackoverflow.com/questions/2100490/floating-point-inaccuracy-examples?rq=1| Floating point inaccuracy examples]</ref> <ref>[https://github.com/gmarkall/PitfallsFP Catastrophic Cancellation: The Pitfalls of Floating Point Arithmetic by Graham Markall]</ref><ref>[https://math.stackexchange.com/questions/2453939/is-this-characteristic-of-tent-map-usually-observed math.stackexchange question: is-this-characteristic-of-tent-map-usually-observed]</ref>
*** porównywanie <ref>[http://stackoverflow.com/questions/10334688/how-dangerous-is-it-to-compare-floating-point-values?rq=1|How dangerous is it to compare floating point values?]</ref>
*** odejmowanie podobnych liczb<ref>[http://stackoverflow.com/questions/28474339/is-it-possible-to-get-0-by-subtracting-two-unequal-floating-point-numbers?rq=1|Is it possible to get 0 by subtracting two unequal floating point numbers?]</ref><ref>[http://fdang.fr/articles/numerical/solve-quadratic-equation.html How to solve quadratic equations numerically by FLORIAN DANG]</ref>
*** błędy zaokrąglania <ref>[http://www.exploringbinary.com/double-rounding-errors-in-floating-point-conversions/ Double Rounding Errors in Floating-Point Conversions By Rick Regan (Published August 11th, 2010)]</ref>
*** przekroczenie limitu ( underflow, overflow)<ref>[http://stackoverflow.com/questions/42277132/when-does-underflow-occur#comment71710792_42277132 stackoverflow question when-does-underflow-occur]</ref><ref>[http://stackoverflow.com/questions/27439503/how-to-define-underflow-for-an-implementationieee754-which-support-subnormal-n stackoverflow question: how-to-define-underflow-for-an-implementationieee754-which-support-subnormal-n]</ref>
*** mnożenie
Efekt:
* zachowanie niezdefiniowane : ( ang. undefined behavior = UB)<ref>[[:w:en:Undefined_behavior|Undefined_behavior w ang. wikipedii]]</ref> <ref>[https://www.nayuki.io/page/undefined-behavior-in-c-and-cplusplus-programs undefined-behavior-in-c-and-cplusplus-programs by Nayuki]</ref>
====Mnożenie ====
<syntaxhighlight lang=c>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
/*
https://math.stackexchange.com/questions/2453939/is-this-characteristic-of-tent-map-usually-observed
*/
/* ------------ constans ---------------------------- */
double m = 2.0; /* parameter of tent map */
double a = 1.0; /* upper bound for randum number generator */
int iMax = 100;
/* ------------------- functions --------------------------- */
/*
tent map
https://en.wikipedia.org/wiki/Tent_map
*/
double f(double x0, double m){
double x1;
if (x0 < 0.5)
x1 = m*x0;
else x1 = m*(1.0 - x0);
return x1;
}
/* random double from 0.0 to a
https://stackoverflow.com/questions/13408990/how-to-generate-random-float-number-in-c
*/
double GiveRandom(double a){
srand((unsigned int)time(NULL));
return (((double)rand()/(double)(RAND_MAX)) * a);
}
int main(void){
int i = 0;
double x = GiveRandom(a); /* x0 = random */
for (i = 0; i<iMax; i++){
printf("i = %3d \t x = %.16f\n",i, x);
x = f(x,m); /* iteration of the tent map */
}
return 0;
}
</syntaxhighlight>
Kompilacja i uruchomienie:
gcc t.c -Wall
./a.out
Wynik:
<pre>
i = 0 x = 0.1720333817284710
i = 1 x = 0.3440667634569419
i = 2 x = 0.6881335269138839
i = 3 x = 0.6237329461722323
i = 4 x = 0.7525341076555354
i = 5 x = 0.4949317846889292
i = 6 x = 0.9898635693778584
i = 7 x = 0.0202728612442833
i = 8 x = 0.0405457224885666
i = 9 x = 0.0810914449771332
i = 10 x = 0.1621828899542663
i = 11 x = 0.3243657799085327
i = 12 x = 0.6487315598170653
i = 13 x = 0.7025368803658694
i = 14 x = 0.5949262392682613
i = 15 x = 0.8101475214634775
i = 16 x = 0.3797049570730451
i = 17 x = 0.7594099141460902
i = 18 x = 0.4811801717078197
i = 19 x = 0.9623603434156394
i = 20 x = 0.0752793131687213
i = 21 x = 0.1505586263374425
i = 22 x = 0.3011172526748851
i = 23 x = 0.6022345053497702
i = 24 x = 0.7955309893004596
i = 25 x = 0.4089380213990808
i = 26 x = 0.8178760427981615
i = 27 x = 0.3642479144036770
i = 28 x = 0.7284958288073540
i = 29 x = 0.5430083423852921
i = 30 x = 0.9139833152294159
i = 31 x = 0.1720333695411682
i = 32 x = 0.3440667390823364
i = 33 x = 0.6881334781646729
i = 34 x = 0.6237330436706543
i = 35 x = 0.7525339126586914
i = 36 x = 0.4949321746826172
i = 37 x = 0.9898643493652344
i = 38 x = 0.0202713012695312
i = 39 x = 0.0405426025390625
i = 40 x = 0.0810852050781250
i = 41 x = 0.1621704101562500
i = 42 x = 0.3243408203125000
i = 43 x = 0.6486816406250000
i = 44 x = 0.7026367187500000
i = 45 x = 0.5947265625000000
i = 46 x = 0.8105468750000000
i = 47 x = 0.3789062500000000
i = 48 x = 0.7578125000000000
i = 49 x = 0.4843750000000000
i = 50 x = 0.9687500000000000
i = 51 x = 0.0625000000000000
i = 52 x = 0.1250000000000000
i = 53 x = 0.2500000000000000
i = 54 x = 0.5000000000000000
i = 55 x = 1.0000000000000000
i = 56 x = 0.0000000000000000
i = 57 x = 0.0000000000000000
i = 58 x = 0.0000000000000000
i = 59 x = 0.0000000000000000
i = 60 x = 0.0000000000000000
i = 61 x = 0.0000000000000000
i = 62 x = 0.0000000000000000
i = 63 x = 0.0000000000000000
i = 64 x = 0.0000000000000000
i = 65 x = 0.0000000000000000
i = 66 x = 0.0000000000000000
i = 67 x = 0.0000000000000000
i = 68 x = 0.0000000000000000
i = 69 x = 0.0000000000000000
i = 70 x = 0.0000000000000000
i = 71 x = 0.0000000000000000
i = 72 x = 0.0000000000000000
i = 73 x = 0.0000000000000000
i = 74 x = 0.0000000000000000
i = 75 x = 0.0000000000000000
i = 76 x = 0.0000000000000000
i = 77 x = 0.0000000000000000
i = 78 x = 0.0000000000000000
i = 79 x = 0.0000000000000000
i = 80 x = 0.0000000000000000
i = 81 x = 0.0000000000000000
i = 82 x = 0.0000000000000000
i = 83 x = 0.0000000000000000
i = 84 x = 0.0000000000000000
i = 85 x = 0.0000000000000000
i = 86 x = 0.0000000000000000
i = 87 x = 0.0000000000000000
i = 88 x = 0.0000000000000000
i = 89 x = 0.0000000000000000
i = 90 x = 0.0000000000000000
i = 91 x = 0.0000000000000000
i = 92 x = 0.0000000000000000
i = 93 x = 0.0000000000000000
i = 94 x = 0.0000000000000000
i = 95 x = 0.0000000000000000
i = 96 x = 0.0000000000000000
i = 97 x = 0.0000000000000000
i = 98 x = 0.0000000000000000
i = 99 x = 0.0000000000000000
</pre>
====Porównywanie====
Sprawdźmy czy liczba x jest równa zero :
if (x==0.0)
Czy takie porównanie jest bezpieczne dla liczb zmiennoprzecinkowych ?<ref>[https://stackoverflow.com/questions/20362304/c-floating-point-zero-comparison stackoverflow question: c-floating-point-zero-comparison]</ref>
<syntaxhighlight lang=c>
// gcc c1.c -Wall -lm
#include <math.h> /* isnormal */
#include <float.h>//DBL_MIN
#include <stdio.h>
int main ()
{
double x = 1.0;
int i;
for ( i=0; i < 334; i++)
{
x/=10;
printf ("i = %3d ; x= %.16lf = %e so ", i, x,x);
//
if (x<DBL_MIN) printf ("x < DBL_MIN and ");
else printf ("x > DBL_MIN and ");
//
if (isnormal(x)) printf ("x is normal and ");
else printf ("x is subnormal and ");
//
if (x==0.0) printf ("equal to 0.0\n");
else printf ("not equal to 0.0\n");
}
return 0;
}
</syntaxhighlight>
Wynik :
<pre>
i = 0 ; x= 0.1000000000000000 = 1.000000e-01 so x > DBL_MIN and x is normal and not equal to 0.0
i = 1 ; x= 0.0100000000000000 = 1.000000e-02 so x > DBL_MIN and x is normal and not equal to 0.0
i = 2 ; x= 0.0010000000000000 = 1.000000e-03 so x > DBL_MIN and x is normal and not equal to 0.0
i = 3 ; x= 0.0001000000000000 = 1.000000e-04 so x > DBL_MIN and x is normal and not equal to 0.0
i = 4 ; x= 0.0000100000000000 = 1.000000e-05 so x > DBL_MIN and x is normal and not equal to 0.0
i = 5 ; x= 0.0000010000000000 = 1.000000e-06 so x > DBL_MIN and x is normal and not equal to 0.0
i = 6 ; x= 0.0000001000000000 = 1.000000e-07 so x > DBL_MIN and x is normal and not equal to 0.0
i = 7 ; x= 0.0000000100000000 = 1.000000e-08 so x > DBL_MIN and x is normal and not equal to 0.0
i = 8 ; x= 0.0000000010000000 = 1.000000e-09 so x > DBL_MIN and x is normal and not equal to 0.0
i = 9 ; x= 0.0000000001000000 = 1.000000e-10 so x > DBL_MIN and x is normal and not equal to 0.0
i = 10 ; x= 0.0000000000100000 = 1.000000e-11 so x > DBL_MIN and x is normal and not equal to 0.0
i = 11 ; x= 0.0000000000010000 = 1.000000e-12 so x > DBL_MIN and x is normal and not equal to 0.0
i = 12 ; x= 0.0000000000001000 = 1.000000e-13 so x > DBL_MIN and x is normal and not equal to 0.0
i = 13 ; x= 0.0000000000000100 = 1.000000e-14 so x > DBL_MIN and x is normal and not equal to 0.0
i = 14 ; x= 0.0000000000000010 = 1.000000e-15 so x > DBL_MIN and x is normal and not equal to 0.0
i = 15 ; x= 0.0000000000000001 = 1.000000e-16 so x > DBL_MIN and x is normal and not equal to 0.0
i = 16 ; x= 0.0000000000000000 = 1.000000e-17 so x > DBL_MIN and x is normal and not equal to 0.0
i = 17 ; x= 0.0000000000000000 = 1.000000e-18 so x > DBL_MIN and x is normal and not equal to 0.0
i = 18 ; x= 0.0000000000000000 = 1.000000e-19 so x > DBL_MIN and x is normal and not equal to 0.0
i = 19 ; x= 0.0000000000000000 = 1.000000e-20 so x > DBL_MIN and x is normal and not equal to 0.0
i = 20 ; x= 0.0000000000000000 = 1.000000e-21 so x > DBL_MIN and x is normal and not equal to 0.0
...
i = 290 ; x= 0.0000000000000000 = 1.000000e-291 so x > DBL_MIN and x is normal and not equal to 0.0
i = 291 ; x= 0.0000000000000000 = 1.000000e-292 so x > DBL_MIN and x is normal and not equal to 0.0
i = 292 ; x= 0.0000000000000000 = 1.000000e-293 so x > DBL_MIN and x is normal and not equal to 0.0
i = 293 ; x= 0.0000000000000000 = 1.000000e-294 so x > DBL_MIN and x is normal and not equal to 0.0
i = 294 ; x= 0.0000000000000000 = 1.000000e-295 so x > DBL_MIN and x is normal and not equal to 0.0
i = 295 ; x= 0.0000000000000000 = 1.000000e-296 so x > DBL_MIN and x is normal and not equal to 0.0
i = 296 ; x= 0.0000000000000000 = 1.000000e-297 so x > DBL_MIN and x is normal and not equal to 0.0
i = 297 ; x= 0.0000000000000000 = 1.000000e-298 so x > DBL_MIN and x is normal and not equal to 0.0
i = 298 ; x= 0.0000000000000000 = 1.000000e-299 so x > DBL_MIN and x is normal and not equal to 0.0
i = 299 ; x= 0.0000000000000000 = 1.000000e-300 so x > DBL_MIN and x is normal and not equal to 0.0
i = 300 ; x= 0.0000000000000000 = 1.000000e-301 so x > DBL_MIN and x is normal and not equal to 0.0
i = 301 ; x= 0.0000000000000000 = 1.000000e-302 so x > DBL_MIN and x is normal and not equal to 0.0
i = 302 ; x= 0.0000000000000000 = 1.000000e-303 so x > DBL_MIN and x is normal and not equal to 0.0
i = 303 ; x= 0.0000000000000000 = 1.000000e-304 so x > DBL_MIN and x is normal and not equal to 0.0
i = 304 ; x= 0.0000000000000000 = 1.000000e-305 so x > DBL_MIN and x is normal and not equal to 0.0
i = 305 ; x= 0.0000000000000000 = 1.000000e-306 so x > DBL_MIN and x is normal and not equal to 0.0
i = 306 ; x= 0.0000000000000000 = 1.000000e-307 so x > DBL_MIN and x is normal and not equal to 0.0
i = 307 ; x= 0.0000000000000000 = 1.000000e-308 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 308 ; x= 0.0000000000000000 = 1.000000e-309 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 309 ; x= 0.0000000000000000 = 1.000000e-310 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 310 ; x= 0.0000000000000000 = 1.000000e-311 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 311 ; x= 0.0000000000000000 = 1.000000e-312 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 312 ; x= 0.0000000000000000 = 1.000000e-313 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 313 ; x= 0.0000000000000000 = 1.000000e-314 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 314 ; x= 0.0000000000000000 = 1.000000e-315 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 315 ; x= 0.0000000000000000 = 1.000000e-316 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 316 ; x= 0.0000000000000000 = 9.999997e-318 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 317 ; x= 0.0000000000000000 = 9.999987e-319 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 318 ; x= 0.0000000000000000 = 9.999889e-320 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 319 ; x= 0.0000000000000000 = 9.999889e-321 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 320 ; x= 0.0000000000000000 = 9.980126e-322 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 321 ; x= 0.0000000000000000 = 9.881313e-323 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 322 ; x= 0.0000000000000000 = 9.881313e-324 so x < DBL_MIN and x is subnormal and not equal to 0.0
i = 323 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 324 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 325 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 326 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 327 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 328 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 329 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 330 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 331 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 332 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
i = 333 ; x= 0.0000000000000000 = 0.000000e+00 so x < DBL_MIN and x is subnormal and equal to 0.0
</pre>
Jak powinno się porównywać liczby zmienno przecinkowe ?<ref>[https://bitbashing.io/comparing-floats.html Comparing Floating-Point Numbers Is Tricky by Matt Kline]</ref>
* wartość bezwględna róznicy : if( abs(a-b) < epsilon) // wrong - don't do this<ref>[https://floating-point-gui.de/errors/comparison/ | floating-point-gui.de : comparison/]</ref>
* if( abs((a-b)/b) < epsilon ) // still not right!
* wartości graniczne
* stałe <ref>[https://stackoverflow.com/questions/9542391/float-double-equality-with-exact-zero stackoverflow question : float-double-equality-with-exact-zero]</ref>
if (fpclassify(x) == FP_ZERO )
lub
if (x == FP_ZERO)
Wartości służące do testo wania porównań :
* wg Michael Borgwardt<ref>[http://floating-point-gui.de/errors/NearlyEqualsTest.java Michael Borgwardt : Nearly Equals Test in java]</ref>
====Sumowanie====
Na ile poważny jest to problem? Spróbujmy przyjrzeć się działaniu, polegającym na 1000-krotnym dodawaniu do liczby wartości 1/3. Oto kod:
<syntaxhighlight lang="c">
#include <stdio.h>
int main ()
{
float a = 0;
int i = 0;
for (;i<1000;i++)
{
a += 1.0/3.0;
}
printf ("%f\n", a);
}
</syntaxhighlight>
Z matematyki wynika, że 1000*(1/3) = 333.(3), podczas gdy komputer wypisze wynik, nieco różniący się od oczekiwanego (w moim przypadku):
<pre>
333.334106
</pre>
Błąd pojawił się na cyfrze części tysięcznej liczby. Nie jest to może poważny błąd, jednak zastanówmy się, czy ten błąd nie będzie się powiększał. Zamieniamy w kodzie ilość iteracji z 1000 na 100 000. Tym razem mój komputer wskazał już nieco inny wynik:
<pre>
33356.554688
</pre>
Błąd przesunął się na cyfrę dziesiątek w liczbie. Tak więc nie należy do końca polegać na prezentacji liczb zmiennoprzecinkowych w komputerze.
====Utrata precyzji====
Utrata precyzji, utrata cyfr znaczących ( ang. Loss of significance, catastrophic cancellation of significant digits)
* sumowanie dużej liczby z małą
* odejmowanie prawie równych liczb<ref>[https://www.nayuki.io/page/numerically-stable-law-of-cosines numerically-stable-law-of-cosines by Nayuki]</ref>
==Biblioteki matematyczne ==
|