C/Zaawansowane operacje matematyczne

Z Wikibooks, biblioteki wolnych podręczników.
< C
Skocz do: nawigacja, szukaj

Prezentacja liczb rzeczywistych w pamięci komputera[edytuj]

W wielu książkach nie ma w ogóle tego tematu. Być może ten temat może wydać Ci się niepotrzebnym, lecz dzięki niemu zrozumiesz, jak komputer radzi sobie z przecinkiem oraz dlaczego niektóre obliczenia dają niezbyt dokładne wyniki.[1]
Na początek trochę teorii. Do przechowywania liczb rzeczywistych przeznaczone są 3 typy: float, double oraz long double. Zajmują one odpowiednio 32, 64 oraz 80 bitów. Wiemy też, że komputer nie ma fizycznej możliwości zapisania przecinka. Spróbujmy teraz zapisać jakąś liczbę wymierną w formie liczb binarnych. Nasza liczba to powiedzmy 4.25. Spróbujmy ją rozbić na sumę potęg dwójki: 4 = 1*22 + 0*21+0*20. Dobra - rozpisaliśmy liczbę 4, ale co z częścią dziesiętną? Skorzystajmy z zasad matematyki - 0.25 = 2-2. Zatem nasza liczba powinna wyglądać tak:

100.01

Ponieważ komputer nie jest w stanie przechować pozycji przecinka, ktoś wpadł na prosty ale sprytny pomysł ustawienia przecinka jak najbliżej początku liczby i tylko mnożenia jej przez odpowiednią potęgę dwójki. Taki sposób przechowywania liczb nazywamy zmiennoprzecinkowym, a proces przekształcania naszej liczby z postaci czytelnej przez człowieka na format zmiennoprzecinkowy nazywamy normalizacją. Wróćmy do naszej liczby - 4.25. W postaci binarnej wygląda ona tak: 100.01, natomiast po normalizacji będzie wyglądała tak: 1.0001*22. W ten sposób w pamięci komputera znajdą się dwie informacje: liczba zakodowana w pamięci z "wirtualnym" przecinkiem oraz numer potęgi dwójki. Te dwie informacje wystarczają do przechowania wartości liczby. Jednak pojawia się inny problem - co się stanie, jeśli np. będziemy chcieli przełożyć liczbę typu \frac{1}{3}? Otóż tutaj wychodzą na wierzch pewne niedociągnięcia komputera w dziedzinie samej matematyki. 1/3 daje w rozwinięciu dziesiętnym 0.(3). Jak zatem zapisać taką liczbę? Otóż nie możemy przechować całego jej rozwinięcia (wynika to z ograniczeń typu danych - ma on niestety skończoną liczbę bitów). Dlatego przechowuje się tylko pewne przybliżenie liczby. Jest ono tym bardziej dokładne im dany typ ma więcej bitów. Zatem do obliczeń wymagających dokładnych danych powinniśmy użyć typu double lub long double. Na szczęście w większości przeciętnych programów tego typu problemy zwykle nie występują. A ponieważ początkujący programista nie odpowiada za tworzenie programów sterujących np. lotem statku kosmicznego, więc drobne przekłamania na odległych miejscach po przecinku nie stanowią większego problemu.

Program wyświetlający wewnętrzną reprezentację liczby podwójnej precyzji:[2]

#include <stdio.h>
 
int main(void)
{
  double a = 1.0 / 3;
  size_t i;
  size_t iMax= sizeof a;
 
  printf("bytes are numbered from 0 to %x\n", (unsigned)iMax -1);
 
  for (i = 0; i < iMax ; ++i) 
  {
    printf("byte number %u  =  %x\n", (unsigned)i, ((unsigned char *)&a)[i]);
  }
 
 
  printf("hex memory representation  of %f is : \n", a);
  for (i = iMax; i>0 ; --i) 
  {
    printf("%x", ((unsigned char *)&a)[i-1]);
  }
  printf(" \n");
  return 0;
}

Daje wynik:

bytes are numbered from 0 to 7
byte number 0  =  55
byte number 1  =  55
byte number 2  =  55
byte number 3  =  55
byte number 4  =  55
byte number 5  =  55
byte number 6  =  d5
byte number 7  =  3f

hex memory representation  of 0.333333 is : 
3fd5555555555555

Obliczenia numeryczne[edytuj]

Obliczenia numeryczne[3] są to obliczenia na liczbach. Ich przeciwieństwem są obliczenia symboliczne wykonywane na symbolach ( zobacz Maxima CAS ) [4]

Uwaga! 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.[5] 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.[6]

Błędy w obliczeniach numerycznych[edytuj]

Typy:[7]

  • błędy zaokrąglania [8]

Przykład[edytuj]

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:

#include <stdio.h>
 
int main ()
{
  float a = 0;
  int i = 0;
  for (;i<1000;i++)
  {
    a += 1.0/3.0;
  }
  printf ("%f\n", a);
}

Z matematyki wynika, że 1000*(1/3) = 333.(3), podczas gdy komputer wypisze wynik, nieco różniący się od oczekiwanego (w moim przypadku):

333.334106

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:

33356.554688

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[edytuj]

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

Epsilon maszynowy[edytuj]

Epsilon maszynowy jest wartością określającą precyzję obliczeń numerycznych wykonywanych na liczbach zmiennoprzecinkowych.[9]

Jest to najmniejsza liczba nieujemna, której dodanie do jedności daje wynik nierówny 1. Innymi słowy, jest to najmniejszy ε, dla którego następujący warunek jest uznawany za niespełniony (przyjmuje wartość fałsz): 1 + ε = 1

Im mniejsza wartość epsilona maszynowego, tym większa jest względna precyzja obliczeń.

Uwaga! Uwaga!
Wartości tej nie należy mylić ze (zwykle dużo niższą) najmniejszą liczbą uznawaną przez maszynę za różną od zera ( np. DBL_MIN dla typu double ).

Obliczmy epsilon dla liczb podwójnej precyzji :

/* 
http://en.wikipedia.org/wiki/Machine_epsilon
The following C program does not actually determine the machine epsilon;
rather, it determines a number within a factor of two (one order of magnitude) 
of the true machine epsilon, using a linear search.
 
gcc m.c -lm -Wall
 
*/
#include <stdio.h>
 
 int main()
 {
    double epsilon = 1.0;
 
    printf( "epsilon;  1 + epsilon\n" );
 
    do 
     {
       printf( "%G\t%.20f\n", epsilon, (1.0 + epsilon) );
       epsilon /= 2.0f;
     }
    // If next epsilon yields 1, then break
    while ((1.0 + (epsilon/2.0)) != 1.0); // 
 
    // because current epsilon is the machine epsilon.
    printf( "\nCalculated Machine epsilon: %G\n", epsilon );
    return 0;
 }

Wynik programu :

epsilon;  1 + epsilon
1               2.00000000000000000000
0.5             1.50000000000000000000
0.25            1.25000000000000000000
0.125           1.12500000000000000000
0.0625          1.06250000000000000000
0.03125         1.03125000000000000000
0.015625        1.01562500000000000000
0.0078125       1.00781250000000000000
0.00390625      1.00390625000000000000
0.00195312      1.00195312500000000000
0.000976562     1.00097656250000000000
0.000488281     1.00048828125000000000
0.000244141     1.00024414062500000000
0.00012207      1.00012207031250000000
6.10352E-05     1.00006103515625000000
3.05176E-05     1.00003051757812500000
1.52588E-05     1.00001525878906250000
7.62939E-06     1.00000762939453125000
3.8147E-06      1.00000381469726562500
1.90735E-06     1.00000190734863281250
9.53674E-07     1.00000095367431640625
4.76837E-07     1.00000047683715820312
2.38419E-07     1.00000023841857910156
1.19209E-07     1.00000011920928955078
5.96046E-08     1.00000005960464477539
2.98023E-08     1.00000002980232238770
1.49012E-08     1.00000001490116119385
7.45058E-09     1.00000000745058059692
3.72529E-09     1.00000000372529029846
1.86265E-09     1.00000000186264514923
9.31323E-10     1.00000000093132257462
4.65661E-10     1.00000000046566128731
2.32831E-10     1.00000000023283064365
1.16415E-10     1.00000000011641532183
5.82077E-11     1.00000000005820766091
2.91038E-11     1.00000000002910383046
1.45519E-11     1.00000000001455191523
7.27596E-12     1.00000000000727595761
3.63798E-12     1.00000000000363797881
1.81899E-12     1.00000000000181898940
9.09495E-13     1.00000000000090949470
4.54747E-13     1.00000000000045474735
2.27374E-13     1.00000000000022737368
1.13687E-13     1.00000000000011368684
5.68434E-14     1.00000000000005684342
2.84217E-14     1.00000000000002842171
1.42109E-14     1.00000000000001421085
7.10543E-15     1.00000000000000710543
3.55271E-15     1.00000000000000355271
1.77636E-15     1.00000000000000177636
8.88178E-16     1.00000000000000088818
4.44089E-16     1.00000000000000044409

Calculated Machine epsilon: 2.22045E-16

real    0m0.005s
user    0m0.000s
sys     0m0.000s


Obliczmy epsilon dla liczb pojedynczej precyzji :

#include <stdio.h>
 
 int main()
 {
    float epsilon = 1.0f;
 
    printf( "epsilon;  1 + epsilon\n" );
 
    do 
     {
       printf( "%G\t%.20f\n", epsilon, (1.0f + epsilon) );
       epsilon /= 2.0f;
     }
    // If next epsilon yields 1, then break
    while ((float)(1.0 + (epsilon/2.0)) != 1.0); // 
 
    // because current epsilon is the machine epsilon.
    printf( "\nCalculated Machine epsilon: %G\n", epsilon );
    return 0;
 }


Wynik :

Calculated Machine epsilon: 1.19209E-07


Obliczmy epsilon dla liczb long double :

#include <stdio.h>
 
 int main()
 {
    long double epsilon = 1.0;
 
    printf( "epsilon;  1 + epsilon\n" );
 
    do 
     {
       printf( "%LG \t %.25LG \n", epsilon, (1.0 + epsilon) );
       epsilon /= 2.0;
     }
    // If next epsilon yields 1, then break
    while ((1.0 + (epsilon/2.0)) != 1.0); // 
 
    // because current epsilon is the machine epsilon.
    printf( "\n Calculated Machine epsilon: %LG\n", epsilon );
    return 0;
 }


Wynik :

 Calculated Machine epsilon: 1.0842E-19

Limity dla obliczeń zmiennoprzecinkowych[edytuj]

Definicje[edytuj]

W pliku float.h są zdefiniowane stałe :

  • DBL_MIN , czyli najmniejszą dodatnia liczbą typu double uznawaną przez maszynę za różną od zera
  • DBL_MAX, czyli największa dodatnia liczbą typu double, która może być używana przez komputer
// gcc -lm -Wall l.c
#include <stdio.h>
#include <math.h> // infinity, nan
#include <float.h>//DBL_MIN
 
int main(void)
{
 
  printf("DBL_MIN = %g \n", DBL_MIN);
  printf("DBL_MAX = %g \n", DBL_MAX);
  printf("INFINITY = %g \n",  INFINITY);
#ifdef NAN 
  printf("NAN= %g \n", NAN );
#endif
  return 0;
}


Wynik działania :

DBL_MIN = 2.22507e-308 
DBL_MAX = 1.79769e+308 
INFINITY = inf 
NAN= nan 

Liczby subnormalne[edytuj]

Korzystając z funkcji isnormal zdefiniowanej w pliku math.h możemy samodzielnie poszukać przybliżenia DBL_MIN i liczby subnormalnej.

/* 
isnormal  example 
ISO C99
http://www.cplusplus.com/reference/cmath/isnormal/
http://www.gnu.org/software/libc/manual/html_node/Floating-Point-Classes.html
http://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
 
compile with: 
gcc -std=c99 s.c -lm
 
run :
./a.out
 
*/
 
#include <stdio.h> /* printf */
#include <math.h> /* isnormal */
 
int TestNumber(double x)
{
  int f; // flag
 
  f= isnormal(x);
  if (f) 
    printf (" = %g ; number is normal \n",x);
    else printf (" = %g ; number is not normal = denormal  \n",x);
 
  return f;
}
 
 
int main()
{
 
  double d ;
  double MinNormal; 
  int flag;
 
  d = 1.0 ; // normal
  flag = TestNumber(d);
  do 
  { 
   MinNormal=d;
   d /=2.0; // look for subnormal 
   flag = TestNumber(d);
 
   }
  while (flag);
 
  printf ("number %f = %g = %e is a approximation of minimal positive double normal \n",MinNormal, MinNormal, MinNormal);
  printf ("number %f = %g = %e is not normal ( subnormal) \n",d, d , d);
 
 return 0;
}

Wynik działania :

number 0.000000 = 2.22507e-308 = 2.225074e-308 is a approximation of minimal positive double normal 
number 0.000000 = 1.11254e-308 = 1.112537e-308 is not normal ( subnormal) 


Biblioteki matematyczne[edytuj]

Standardowa : math.h[edytuj]

Aby móc korzystać z wszystkich dobrodziejstw funkcji matematycznych musimy na początku dołączyć plik math.h:

 #include <math.h>

a w procesie kompilacji (dotyczy kompilatora GCC) musimy dodać flagę "-lm" po nazwie pliku wynikowego[10], czyli na końcu linii :[11]

gcc plik.c -o plik -lm

Funkcje matematyczne, które znajdują się w bibliotece standardowej ( plik libm.a ) możesz znaleźć tutaj. Przy korzystaniu z nich musisz wziąć pod uwagę m.in. to, że biblioteka matematyczna prowadzi kalkulację w oparciu o radiany a nie stopnie.

Stałe matematyczne[edytuj]

W pliku math.h zdefiniowane są pewne stałe, które mogą być przydatne do obliczeń. Są to m.in.:

  • M_E - podstawa logarytmu naturalnego (e, liczba Eulera)
  • M_LOG2E - logarytm o podstawie 2 z liczby e
  • M_LOG10E - logarytm o podstawie 10 z liczby e
  • M_LN2 - logarytm naturalny z liczby 2
  • M_LN10 - logarytm naturalny z liczby 10
  • M_PI - liczba π
  • M_PI_2 - liczba π/2
  • M_PI_4 - liczba π/4
  • M_1_PI - liczba 1/π
  • M_2_PI - liczba 2/π


Możemy to sprawdzić:

grep -i pi /usr/include/math.h

i otrzymamy:

# define M_PI        3.14159265358979323846    /* pi */
# define M_PI_2        1.57079632679489661923    /* pi/2 */
# define M_PI_4        0.78539816339744830962    /* pi/4 */
# define M_1_PI        0.31830988618379067154    /* 1/pi */
# define M_2_PI        0.63661977236758134308    /* 2/pi */
# define M_2_SQRTPI    1.12837916709551257390    /* 2/sqrt(pi) */
# define M_PIl        3.1415926535897932384626433832795029L  /* pi */
# define M_PI_2l    1.5707963267948966192313216916397514L  /* pi/2 */
# define M_PI_4l    0.7853981633974483096156608458198757L  /* pi/4 */
# define M_1_PIl    0.3183098861837906715377675267450287L  /* 1/pi */
# define M_2_PIl    0.6366197723675813430755350534900574L  /* 2/pi */
# define M_2_SQRTPIl    1.1283791670955125738961589031215452L  /* 2/sqrt(pi) */
/* When compiling in strict ISO C compatible mode we must not use the

Liczby zespolone[edytuj]

complex.h[edytuj]

Operacje na liczbach zespolonych są częścią uaktualnionego standardu języka C o nazwie C99, który jest obsługiwany jedynie przez część kompilatorów

Porada Podane tutaj informacje zostały sprawdzone na systemie Gentoo Linux z biblioteką GNU libc w wersji 2.3.5 i kompilatorem GCC w wersji 4.0.2

Dotychczas korzystaliśmy tylko z liczb rzeczywistych, lecz najnowsze standardy języka C umożliwiają korzystanie także z innych liczb - np. z liczb zespolonych.

Aby móc korzystać z liczb zespolonych w naszym programie należy w nagłówku programu umieścić następującą linijkę:

 #include <complex.h>

Wiemy, że liczba zespolona zdeklarowana jest następująco:

z = a+b*i,

gdzie a, b są liczbami rzeczywistymi, a i jest jednostką urojoną

i*i = (-1).

W pliku complex.h liczba i zdefiniowana jest jako I. Zatem wypróbujmy możliwości liczb zespolonych:

 #include <math.h>
 #include <complex.h>
 #include <stdio.h>
 
 int main ()
 {
   float _Complex z = 4+2.5*I;
   printf ("Liczba z: %f+%fi\n", creal(z), cimag (z));
   return 0;
 }

następnie kompilujemy nasz program:

gcc plik1.c -o plik1 -lm

Po wykonaniu naszego programu powinniśmy otrzymać:

Liczba z: 4.00+2.50i

W programie zamieszczonym powyżej użyliśmy dwóch funkcji - creal i cimag.

  • creal - zwraca część rzeczywistą liczby zespolonej
  • cimag - zwraca część urojoną liczby zespolonej


MPC[edytuj]

Opis MPC

Dodatkowe[edytuj]

Źródła[edytuj]

  1. Metody numeryczne - autorzy : Piotr Krzyżanowski i Leszek Plaskota — Uniwersytet Warszawski, Wydział Matematyki, Informatyki i Mechaniki
  2. How to get memory representation-double
  3. numeryczne - Wydziału MIM UW
  4. i ból obliczeń numerycznych -Piotr Krzyżanowski
  5. Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg
  6. Two disasters caused by computer arithmetic errors
  7. [www.sci.brooklyn.cuny.edu/~mate/nml/numanal.pdf|INTRODUCTION TO NUMERICAL ANALYSIS WITH C PROGRAMS by Attila Mate]
  8. [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)]
  9. Praktyczne wyznaczanie precyzji arytmetyki - autorzy : Piotr Krzyżanowski i Leszek Plaskota — Uniwersytet Warszawski, Wydział Matematyki, Informatyki i Mechaniki
  10. Stackoverflow : Undefined reference to `pow' and `floor'
  11. I'm using math.h and the library link -lm, but “undefined reference to `pow'” still happening