C/Zaawansowane operacje matematyczne

< C

Biblioteki

edytuj
  • Biblioteka math.h
  • complex.h - liczby zespolone
  • liczby rzeczywiste
    • float.h - liczny zmiennoprzecinkowe
    • fenv.h - środowiski zmiennoprzecinkowe
  • liczby całkowite
    • limits.h - charakterystyka typów liczb calkowitych
    • inittypes.h - konwersja liczb całkowitych
    • stdckdiniot.h - kontrola arytmetyki liczb całkowitych
  • atomic.h - typy atomowe
  • stdbit.h - marzędzia bitowe i bajtowe
  • stdbool.h - wartości i typy logiczne
  • tgmath.h - funkcje matematyczne niezależne od typu


Liczby

edytuj
Liczby 0–F16 w systemie
dziesiętnym, ósemkowym i dwójkowym
  0hex   =   0dec   =   0oct     0     0     0     0  
1hex = 1dec = 1oct 0 0 0 1
2hex = 2dec = 2oct 0 0 1 0
3hex = 3dec = 3oct 0 0 1 1
4hex = 4dec = 4oct 0 1 0 0
5hex = 5dec = 5oct 0 1 0 1
6hex = 6dec = 6oct 0 1 1 0
7hex = 7dec = 7oct 0 1 1 1
8hex = 8dec = 10oct 1 0 0 0
9hex = 9dec = 11oct 1 0 0 1
Ahex = 10dec = 12oct 1 0 1 0
Bhex = 11dec = 13oct 1 0 1 1
Chex = 12dec = 14oct 1 1 0 0
Dhex = 13dec = 15oct 1 1 0 1
Ehex = 14dec = 16oct 1 1 1 0
Fhex = 15dec = 17oct 1 1 1 1

Rodzaje liczb

edytuj

Dodatkowo rodzaj liczby definiują :


Jak widać nie ma :

  • Liczb niewymiernych ( chyba że korzystamy z obliczeń symbolicznych, np. program Maxima CAS )

liczby binarne

edytuj

Sposoby :


Binarna stała całkowita z użyciem rozszerzenia gcc [4]

  • prefix ‘0b’ lub ‘0B’
  • sekwencja cyfr ‘0’ lub ‘1’
  • suffix ‘L’ lub ‘UL’

Następujące zapisy są identyczne:

     i =       42; // decimal
     i =     0x2a; // hexadecimal
     i =      052; // octal
     i = 0b101010; // binary


int i = 1 << 9; /* binary 1 followed by 9 binary 0's */ [5]


Liczby szestnastkowe ( hexadecymalne)

edytuj

Liczba szesnastkowa:

  • jest poprzedzana przez przedrostek „0x” lub „0X”
  • jest wartością całkowitą i można ją przechowywać w typie całkowitym typów danych (char, short lub int)


Przykłady:

  • „0x64” odpowiada 100 w systemie dziesiętnym
  • wartość szesnastkowa 3E8 jest reprezentowana w c jako 0x3E8

Każdy bajt ( 8 bitów) to 2-cyfrowa liczba szesnastkowa


Maksymalna liczba całkowita bez znaku zajmująca n bajtów:

  • 1 bajt = 8 bitów:  
  • 2 bajty = 16 bitów  
  • 4 bajty = 32 bity  
  • 8 bajtów = 64 bity  

Prezentacja liczb rzeczywistych w pamięci komputera

edytuj

Liczby rzeczywisty

edytuj
 
Format liczby podwójnej precyzji

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 [6]oraz dlaczego niektóre obliczenia dają niezbyt dokładne wyniki.[7]
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  ? 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:[8][9]

#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


0 01111111101 01010101010101010101010101010101010101010101010101012 = 3FD5 5555 5555 555516 ≙ +2−2 × (1 + 2−2 + 2−4 + ... + 2−52) ≈ 1/3


Given the hexadecimal representation 3FD5 5555 5555 555516,
  Sign = 0
  Exponent = 3FD16 = 1021
  Exponent Bias = 1023 (constant value; see above)
  Fraction = 5 5555 5555 555516
  Value = 2(Exponent − Exponent Bias) × 1.Fraction – Note that Fraction must not be converted to decimal here
        = 2−2 × (15 5555 5555 555516 × 2−52)
        = 2−54 × 15 5555 5555 555516
        = 0.333333333333333314829616256247390992939472198486328125
        ≈ 1/3



Zobacz również:

Liczby całkowite

edytuj

Liczba całkowitej bez znaku (ang. unsigned char lub uint8_t ) zajmującej 1 bajt, czyli 8 bitów.


Ta tabela ilustruje wartość

  

Oznaczenia

  • MSB = najbardziej znaczący bit
  • LSB = najmniej znaczący bit
 
Liczba dziesiętna bez znaku 149 w postaci binarnej z zaznaczonym LSB
 
Liczba dziesiętna bez znaku 149 w postaci binarnej z zaznaczonym MSB

 

Bit position i 7 6 5 4 3 2 1 0
Binary digit   1 0 0 1 0 1 0 1
Bit weight = ( 2i ) 27 26 25 24 23 22 21 20
Bit position label MSB LSB

Liczba całkowita bez znaku zajmująca 2 bajty = 16 bitów ( uint16_t )


Wartość

  • dziesiętna (decymalna) = 769
  • dwójkowa (binarna) = 1100000001
  • hexadecymalna = 0x0301

 


 +---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+
 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |  
 +---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+


Zadania:

  • konwersja dwóch wartości 8-bitowych na jedną wartość 16-bitową w C[11]
  • dzielenie całkowitej liczby 16 bitowej bez znaku na 2 bajty[12]

Konwersje liczb

edytuj
/*

gcc b.c -Wall -Wextra


*/


#include <stdio.h>
#include <stdlib.h>


/*
https://stackoverflow.com/questions/5488377/converting-an-integer-to-binary-in-c


If you want to transform a number into another number (not number to string of characters), and you can do with a small range (0 to 1023 for implementations with 32-bit integers), you don't need to add char* to the solution. pmg

*/
unsigned unsigned_int_to_bin(unsigned int k) {
    if (k == 0) return 0;
    if (k == 1) return 1;                       /* optional */
    return (k % 2) + 10 * unsigned_int_to_bin(k / 2);
}



int main(void){


	unsigned int u = 1019; 
	if (u>1023) {fprintf(stderr, "input number decimal u = %d is too big; close the program ! \n", u); return 1;}
	
	unsigned int b = unsigned_int_to_bin(u);
	fprintf(stdout, "decimal u = %d as a binary unsigned number = %d\n", u, b);
	



	return 0;

}

Wynik:

gcc b.c -Wall -Wextra
./a.out
decimal u = 1019 as a binary unsigned number = 1111111011



Konwersja zapisu matematycznego na program komputerowy w języku C

edytuj

sumowanie

edytuj

W zapisie matematycznym do przedstawiania w zwarty sposób sumowania wielu podobnych wyrazów ( serii) stosuje się symbol sumowania  , wywodzący się z wielkiej greckiej litery sigma. Jest on zdefiniowany jako:

 

gdzie

  •   oznacza indeks sumowania
  •   to zmienna przedstawiająca każdy kolejny wyraz w szeregu
  •   to dolna granica sumowania
  •   to górna granica sumowania.

Wyrażenie „i = m” pod symbolem sumowania oznacza, że indeks   rozpoczyna się od wartości   Następnie dla każdego kolejnego wyrazu indeks   jest zwiększany o 1, aż   osiągnie wartość   (tj.  ), który jest ostatnim wyrazem sumowania.


#include <stdio.h>

int summation(const int m,  const int n )
{
   int s = 0;
   for(int i=m; i<=n; ++i)
   {
     s += i;
   }
   return s;
}



int main()
{
   int sum;
   int m = 1; // lower index
   int n = 100; // upper index
   
   sum = summation(m,n);
   
   printf("sum of integer numbers from %d to %d is %d\n", m, n, sum);
   return 0; 
 }
 gcc s.c -Wall -Wextra
 ./a.out
 sum of integer numbers from 1 to 100 is 5050


Suma sum:

 

Sumy są obliczane od prawej do lewej. Oznacza to, że sigma po lewej stronie jest najwyżej zagnieżdżoną pętlą: [13]

    
sum = 0;
for (i =1; i<=4; i++)
    for (j=1; j <=2; j++)
         sum += i * j* j;

produkt ( iloczyn)

edytuj

W zapisie matematycznym do przedstawiania w zwarty sposób iloczynu wielu podobnych wyrazów ( serii) stosuje się symbol iloczynu  , wywodzący się z wielkiej greckiej litery pi. Jest on zdefiniowany jako:

 


#include <stdio.h>

int product(const int m,  const int n )
{
   int p = 1;
   for(int i=m; i<=n; ++i)
   {
     p *= i;
   }
   return p;
}



int main()
{
   int p;
   int m = 1; // lower index
   int n = 10; // upper index
   
    p = product(m,n);
   
   printf("product of integer numbers from %d to %d is %d\n", m, n, p);
   
}

Dla zakresu [1,10] program działa poprawnie

gcc p.c -Wall -Wextra
./a.out
product of integer numbers from 1 to 10 is 3628800

Dla zakresu [1,100] program daje dziwny wynik:

gcc p.c -Wall -Wextra
./a.out
product of integer numbers from 1 to 100 is 0


Wartość iloczynu możemy obliczyć za pomoca wolfram alfa.

Product[k, {k, 1, 100}] = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Liczba ma 159 cyfr, to więcej niż zakres dostępnych standardowych liczb całkowitych. Przekroczenie zakresu powoduje błędny wynik.

Jeśli dolny zakres wynosi 1 to wynik jest równy silni ( ang factorial). Wartość silnii rośnie szybciej niż wzrost funkcji wykładniczej.

 


Sprawdzenie czy nie będzie przekroczenia zakresu liczb powinniśmy wykonać przed wykonaniem działania:

#include <stdio.h>
#include <limits.h> // INT_MAX
/*
http://c-faq.com/misc/sd26.html
functions for ``careful'' multiplication.

*/

int chkmul(int a, int b)
{
	int sign = 1;
	if(a == 0 || b == 0) return 0;
	if(a < 0) { a = -a; sign = -sign; }
	if(b < 0) { b = -b; sign = -sign; }
	if(INT_MAX / b < a) {
		fprintf(stdout,"int overflow\t");
		return (sign > 0) ? INT_MAX : INT_MIN;
	}
	
	
	return sign*a*b;
}


int product(const int m,  const int n )
{
   int p = 1;
   for(int i=m; i<=n; ++i)
   {
     fprintf(stdout,"i = %d \t", i);
     p = chkmul(p,i);
     fprintf(stdout,"p = %d \n", p);
     
   }
   return p;
}



int main()
{
   int p;
   int m = 1; // lower index
   int n = 100; // upper index
   
    p = product(m,n);
   
   printf("product of integer numbers from %d to %d is %d\n", m, n, p);
   
}


gcc p.c -Wall -Wextra
./a.out
i = 1 	p = 1 
i = 2 	p = 2 
i = 3 	p = 6 
i = 4 	p = 24 
i = 5 	p = 120 
i = 6 	p = 720 
i = 7 	p = 5040 
i = 8 	p = 40320 
i = 9 	p = 362880 
i = 10 	p = 3628800 
i = 11 	p = 39916800 
i = 12 	p = 479001600 
i = 13 	int overflow	p = 2147483647 
i = 14 	int overflow	p = 2147483647 
i = 15 	int overflow	p = 2147483647 
i = 16 	int overflow	p = 2147483647 
i = 17 	int overflow	p = 2147483647 
i = 18 	int overflow	p = 2147483647 
i = 19 	int overflow	p = 2147483647 
i = 20 	int overflow	p = 2147483647 
i = 21 	int overflow	p = 2147483647 
i = 22 	int overflow	p = 2147483647 
i = 23 	int overflow	p = 2147483647 
i = 24 	int overflow	p = 2147483647 
i = 25 	int overflow	p = 2147483647 
i = 26 	int overflow	p = 2147483647 
i = 27 	int overflow	p = 2147483647 
i = 28 	int overflow	p = 2147483647 
i = 29 	int overflow	p = 2147483647 
i = 30 	int overflow	p = 2147483647 
i = 31 	int overflow	p = 2147483647 
i = 32 	int overflow	p = 2147483647 
i = 33 	int overflow	p = 2147483647 
i = 34 	int overflow	p = 2147483647 
i = 35 	int overflow	p = 2147483647 
i = 36 	int overflow	p = 2147483647 
i = 37 	int overflow	p = 2147483647 
i = 38 	int overflow	p = 2147483647 
i = 39 	int overflow	p = 2147483647 
i = 40 	int overflow	p = 2147483647 
i = 41 	int overflow	p = 2147483647 
i = 42 	int overflow	p = 2147483647 
i = 43 	int overflow	p = 2147483647 
i = 44 	int overflow	p = 2147483647 
i = 45 	int overflow	p = 2147483647 
i = 46 	int overflow	p = 2147483647 
i = 47 	int overflow	p = 2147483647 
i = 48 	int overflow	p = 2147483647 
i = 49 	int overflow	p = 2147483647 
i = 50 	int overflow	p = 2147483647 
i = 51 	int overflow	p = 2147483647 
i = 52 	int overflow	p = 2147483647 
i = 53 	int overflow	p = 2147483647 
i = 54 	int overflow	p = 2147483647 
i = 55 	int overflow	p = 2147483647 
i = 56 	int overflow	p = 2147483647 
i = 57 	int overflow	p = 2147483647 
i = 58 	int overflow	p = 2147483647 
i = 59 	int overflow	p = 2147483647 
i = 60 	int overflow	p = 2147483647 
i = 61 	int overflow	p = 2147483647 
i = 62 	int overflow	p = 2147483647 
i = 63 	int overflow	p = 2147483647 
i = 64 	int overflow	p = 2147483647 
i = 65 	int overflow	p = 2147483647 
i = 66 	int overflow	p = 2147483647 
i = 67 	int overflow	p = 2147483647 
i = 68 	int overflow	p = 2147483647 
i = 69 	int overflow	p = 2147483647 
i = 70 	int overflow	p = 2147483647 
i = 71 	int overflow	p = 2147483647 
i = 72 	int overflow	p = 2147483647 
i = 73 	int overflow	p = 2147483647 
i = 74 	int overflow	p = 2147483647 
i = 75 	int overflow	p = 2147483647 
i = 76 	int overflow	p = 2147483647 
i = 77 	int overflow	p = 2147483647 
i = 78 	int overflow	p = 2147483647 
i = 79 	int overflow	p = 2147483647 
i = 80 	int overflow	p = 2147483647 
i = 81 	int overflow	p = 2147483647 
i = 82 	int overflow	p = 2147483647 
i = 83 	int overflow	p = 2147483647 
i = 84 	int overflow	p = 2147483647 
i = 85 	int overflow	p = 2147483647 
i = 86 	int overflow	p = 2147483647 
i = 87 	int overflow	p = 2147483647 
i = 88 	int overflow	p = 2147483647 
i = 89 	int overflow	p = 2147483647 
i = 90 	int overflow	p = 2147483647 
i = 91 	int overflow	p = 2147483647 
i = 92 	int overflow	p = 2147483647 
i = 93 	int overflow	p = 2147483647 
i = 94 	int overflow	p = 2147483647 
i = 95 	int overflow	p = 2147483647 
i = 96 	int overflow	p = 2147483647 
i = 97 	int overflow	p = 2147483647 
i = 98 	int overflow	p = 2147483647 
i = 99 	int overflow	p = 2147483647 
i = 100 	int overflow	p = 2147483647 
product of integer numbers from 1 to 100 is 2147483647

Widać że każdy produkt powyżej górnego zakresu=12 jest błędny. Rozwiązaniem jest:

  • użycie bibliotek o dowolnej precyzji ( GMP )
  • użycie typu zmiennoprzecinkowego (double )
#include <stdio.h>

double double_product(const int m,  const int n )
{
   double p = 1.0;
   for(int i=m; i<=n; ++i)
   {
     p*=i;
     
   }
   return p;
}


int main()
{
   double p;
   int m = 1; // lower index
   int n = 100; // upper index
   
    p = double_product(m,n);
   
   printf("product of integer numbers from %d to %d is %.16e\n", m, n, p);
   
}
gcc d.c -Wall -Wextra
./a.out
product of integer numbers from 1 to 100 is 9.3326215443944102e+157

Jak widać pierwsze 17 cyfr znaczących się zgadza. Błąd wynosi około 10^140 czyli jest mniejszy niż 1%

zaawansowane algorytmy

edytuj
  • Jak tworzyć matematyczne oprogramowanie w C?[14]
  • drukowanie liczb algebraicznych[15]

Obliczenia numeryczne

edytuj

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


Epsilon maszynowy

edytuj

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

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

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.

---

The difference between 1 and the least value greater than 1 that is representable in the given floating-point type, b1-p.
-------------------------------
http://stackoverflow.com/questions/1566198/how-to-portably-get-dbl-epsilon-in-c-c

gcc m.c -Wall

./a.out
*/


#include <stdio.h>
#include <float.h> // DBL_EPSILON
 
 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 calculated machine epsilon.
    printf( "\nCalculated Machine epsilon: %G\n", epsilon );
    
    
    //check value from float.h , Steve Jessop
    if ((1.0 + DBL_EPSILON) != 1.0 
        && 
       (1.0 + DBL_EPSILON/2) == 1.0)
     printf("DBL_EPSILON = %g \n", DBL_EPSILON);
     else printf("DBL_EPSILON is not good !!! \n");


    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
DBL_EPSILON = 2.22045e-16 


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

Porównaj

wyjątki

edytuj

Wyjatki

  • programowe
  • sprzętowe

Przykłady

  • A floating-point exception (FPE) = Wyjątek zmiennoprzecinkowy nie jest wyjątkiem programowym. Powstaje na poziomie mikroprocesora lub ISA. FPE może spowodować wysłąnie sygnału o nazwie SIGFPE, z którym można sobie poradzić, ale nie za pomocą C[22]


Wyjątki zmiennoprzecinkowe są kontrolowane przez kod biblioteki w C99, a nie przez flagi kompilatora[23]


    
#include <fenv.h>
#include <math.h>
#include <stdio.h>

#define PRINTEXC(ex, val) printf(#ex ": %s\n", (val & ex) ? "set" : "unset");

double foo(double a, double b) { return sin(a) / b; }

int main()
{
    int e;
    double x;

    feclearexcept(FE_ALL_EXCEPT);

    x = foo(1.2, 3.1);

    e = fetestexcept(FE_ALL_EXCEPT);
    PRINTEXC(FE_DIVBYZERO, e);
    PRINTEXC(FE_INEXACT, e);
    PRINTEXC(FE_INVALID, e);
    PRINTEXC(FE_OVERFLOW, e);
    PRINTEXC(FE_UNDERFLOW, e);

    putchar('\n');

    feclearexcept(FE_ALL_EXCEPT);

    x += foo(1.2, 0.0);

    e = fetestexcept(FE_ALL_EXCEPT);
    PRINTEXC(FE_DIVBYZERO, e);
    PRINTEXC(FE_INEXACT, e);
    PRINTEXC(FE_INVALID, e);
    PRINTEXC(FE_OVERFLOW, e);
    PRINTEXC(FE_UNDERFLOW, e);
    return lrint(x);
}

Wynik:

    

FE_DIVBYZERO: unset
FE_INEXACT: set
FE_INVALID: unset
FE_OVERFLOW: unset
FE_UNDERFLOW: unset

FE_DIVBYZERO: set
FE_INEXACT: set
FE_INVALID: unset
FE_OVERFLOW: unset
FE_UNDERFLOW: unset


Jest też możliwe użycie sygnałów ( zobacz fenv.h):

    
#pragma STDC FENV_ACCESS on

#define _GNU_SOURCE
#include <fenv.h>

int main()
{
#ifdef FE_NOMASK_ENV
    fesetenv(FE_NOMASK_ENV);
#endif

    // ...
}

Limity dla obliczeń

edytuj

zmiennoprzecinkowych

edytuj

Definicje W pliku float.h są zdefiniowane stałe :[24]

  • DBL_MIN , czyli najmniejszą dodatnia liczbą typu double uznawaną przez maszynę za różną od zera [25]
  • DBL_MAX, czyli największa dodatnia liczbą typu double, która może być używana przez komputer

W pliku math.h są zdefiniowane :

// 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 


całkowitych

edytuj
/*

gcc l.c -lm -Wall
./a.out


http://stackoverflow.com/questions/29592898/do-long-long-and-long-have-same-range-in-c-in-64-bit-machine
*/
#include <stdio.h>
#include <math.h> // M_PI; needs -lm also 
#include <limits.h> // INT_MAX, http://pubs.opengroup.org/onlinepubs/009695399/basedefs/limits.h.html





int main(){

double lMax;


 lMax = log2(INT_MAX);
 printf("INT_MAX \t= %25d ; lMax = log2(INT_MAX) \t= %.0f \n",INT_MAX,  lMax); 

 lMax = log2(UINT_MAX);
 printf("UINT_MAX \t= %25u ; lMax = log2(UINT_MAX) \t= %.0f \n", UINT_MAX,  lMax); 

 lMax = log2(LONG_MAX);
 printf("LONG_MAX \t= %25ld ; lMax = log2(LONG_MAX) \t= %.0f \n",LONG_MAX,  lMax); 


 lMax = log2(ULONG_MAX);
 printf("ULONG_MAX \t= %25lu ; lMax = log2(ULONG_MAX) \t= %.0f \n",ULONG_MAX,  lMax); 

 lMax = log2(LLONG_MAX);
 printf("LLONG_MAX \t= %25lld ; lMax = log2(LLONG_MAX) \t= %.0f \n",LLONG_MAX, lMax); 

 lMax = log2(ULLONG_MAX);
 printf("ULLONG_MAX \t= %25llu ; lMax = log2(ULLONG_MAX) \t= %.0f \n",ULLONG_MAX, lMax); 





return 0;
}

Wynik :

INT_MAX 	=                2147483647 ; lMax = log2(INT_MAX) 	= 31 
UINT_MAX 	=                4294967295 ; lMax = log2(UINT_MAX) 	= 32 
LONG_MAX 	=       9223372036854775807 ; lMax = log2(LONG_MAX) 	= 63 
ULONG_MAX 	=      18446744073709551615 ; lMax = log2(ULONG_MAX) 	= 64 
LLONG_MAX 	=       9223372036854775807 ; lMax = log2(LLONG_MAX) 	= 63 
ULLONG_MAX 	=      18446744073709551615 ; lMax = log2(ULLONG_MAX) 	= 64 


Dla typów o stałej szerokości ( ang. Fixed width integer types ):

/*

gcc h.c -Wall -Wextra
./a.out

printf
------------------------------------ 
X ( specifier character): 
	input  = argument of type ‘unsigned int’
	output =  Unsigned hexadecimal integer (uppercase) without 0X prefix



# (  flags sub-specifier)	Used with X specifiers  : the value is preceeded with 0X respectively for values different than zero


*/
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h> // PRIu32 = all those nifty new format specifiers for the intN_t types and their brethren


int main(void)
{
	
	uint8_t  a = 0XFF; 		// 1 byte  =  8 bits
	uint16_t b = 0XFFFF;  		// 2 bytes = 16 bits
	uint32_t c = 0XFFFFFFFF;	// 4 bytes = 32 bits
	uint64_t d = 0XFFFFFFFFFFFFFFFF;// 8 bytes = 64 bits
	
	int bytes = 1;
	int bits = 8* bytes;
	
	
	printf("uint%d_t \ta ( %d byte  =  %d bits)\t = %#X (hexadecimal)\t\t\t= %d \t\t\t = 2^%d - 1  (decimal) \n", bits, bytes, bits, a, a, bits);
	bytes = 2; bits = 8 * bytes;
	printf("uint%d_t\tb ( %d byte  =  %d bits)\t = %#X (hexadecimal)\t\t\t= %d\t\t\t = 2^%d - 1  (decimal) \n", bits, bytes, bits, b, b, bits);
	bytes = 4; bits = 8 * bytes;
	printf("uint%d_t\tc ( %d byte  =  %d bits)\t = %#X (hexadecimal)\t\t= %"PRIu32" \t\t = 2^%d - 1  (decimal) \n", bits, bytes, bits, c, c, bits);
	bytes = 8; bits = 8 * bytes;
	printf("uint%d_t\tc ( %d byte  =  %d bits)\t = %#lX (hexadecimal)\t= %"PRIu64" \t = 2^%d - 1  (decimal) \n", bits, bytes, bits, d, d, bits);
	
	
	return 0;
}

Wynik:

uint8_t 	a ( 1 byte  =  8 bits)	 = 0XFF (hexadecimal)			= 255 			 = 2^8 - 1  (decimal) 
uint16_t	b ( 2 byte  =  16 bits)	 = 0XFFFF (hexadecimal)			= 65535			 = 2^16 - 1  (decimal) 
uint32_t	c ( 4 byte  =  32 bits)	 = 0XFFFFFFFF (hexadecimal)		= 4294967295 		 = 2^32 - 1  (decimal) 
uint64_t	c ( 8 byte  =  64 bits)	 = 0XFFFFFFFFFFFFFFFF (hexadecimal)	= 18446744073709551615 	 = 2^64 - 1  (decimal) 


Limity dla typów całkowitych o stałej szerokości
Typ Znak Bits Bytes Minimum Maximum
int8_t Signed 8 1 −27 = −128 27 − 1 = 127
uint8_t Unsigned 8 1 0 28 − 1 = 255
int16_t Signed 16 2 −215 = −32,768 215 − 1 = 32,767
uint16_t Unsigned 16 2 0 216 − 1 = 65,535
int32_t Signed 32 4 −231 = −2,147,483,648 231 − 1 = 2,147,483,647
uint32_t Unsigned 32 4 0 232 − 1 = 4,294,967,295
int64_t Signed 64 8 −263 = −9,223,372,036,854,775,808 263 − 1 = 9,223,372,036,854,775,807
uint64_t Unsigned 64 8 0 264 − 1 = 18,446,744,073,709,551,615
Przekroczenie zakresu liczb całkowitych
edytuj

Przekroczenie zakresu liczb całkowitych ( ang. integer overflow ) [26] może dotyczyć liczb całkowitych :[27]

  • bez znaku ( " Unsigned integers are defined to wrap around. " )
  • ze znakiem ( powoduje zachowanie niezdefiniowane - może to powodować Złe Rzeczy czyli zagrożenie bezpieczeństwa komputera [28] )
#include <stdio.h>

/* 

a signed integer overflow is undefined behaviour in C 
check b^i 

to compile : 
  gcc i.c -Wall
to run : 

  ./a.out

*/

int main() {

int i;
int b=2; // base 
int p=1; // p = b^i

for ( i=0 ; i<40; i++){
     
     printf(" b^i = %d ^ %d = %d \n", b, i, p);
     p *= b;
     
}


return 0; 
}


Program kompiluje się i uruchamia bez komunikatów o błędach ale wynik nie jest taki jak naiwnie moglibyśmy się spodziewać :


 b^i = 2 ^ 0 = 1 
 b^i = 2 ^ 1 = 2 
 b^i = 2 ^ 2 = 4 
 b^i = 2 ^ 3 = 8 
 b^i = 2 ^ 4 = 16 
 b^i = 2 ^ 5 = 32 
 b^i = 2 ^ 6 = 64 
 b^i = 2 ^ 7 = 128 
 b^i = 2 ^ 8 = 256 
 b^i = 2 ^ 9 = 512 
 b^i = 2 ^ 10 = 1024 
 b^i = 2 ^ 11 = 2048 
 b^i = 2 ^ 12 = 4096 
 b^i = 2 ^ 13 = 8192 
 b^i = 2 ^ 14 = 16384 
 b^i = 2 ^ 15 = 32768 
 b^i = 2 ^ 16 = 65536 
 b^i = 2 ^ 17 = 131072 
 b^i = 2 ^ 18 = 262144 
 b^i = 2 ^ 19 = 524288 
 b^i = 2 ^ 20 = 1048576 
 b^i = 2 ^ 21 = 2097152 
 b^i = 2 ^ 22 = 4194304 
 b^i = 2 ^ 23 = 8388608 
 b^i = 2 ^ 24 = 16777216 
 b^i = 2 ^ 25 = 33554432 
 b^i = 2 ^ 26 = 67108864 
 b^i = 2 ^ 27 = 134217728 
 b^i = 2 ^ 28 = 268435456 
 b^i = 2 ^ 29 = 536870912 
 b^i = 2 ^ 30 = 1073741824 
 b^i = 2 ^ 31 = -2147483648 
 b^i = 2 ^ 32 = 0 
 b^i = 2 ^ 33 = 0 
 b^i = 2 ^ 34 = 0 
 b^i = 2 ^ 35 = 0 
 b^i = 2 ^ 36 = 0 
 b^i = 2 ^ 37 = 0 
 b^i = 2 ^ 38 = 0 
 b^i = 2 ^ 39 = 0 

Na podstawie wyniku możemy ocenić że zmienna int jest typu 32 bitowego , ponieważ obliczenia do 2^30 są poprawne.

Dla liczb bez znaku przekroczenie zakresu powoduje inny efekt ( modulo ) :

#include <stdio.h>

/* 

 
 Unsigned integers are defined to wrap around.
 "When you work with unsigned types, modular arithmetic (also known as "wrap around" behavior) is taking place."
http://stackoverflow.com/questions/7221409/is-unsigned-integer-subtraction-defined-behavior

*/


int main() {

unsigned int i;
unsigned int old=0; // 
unsigned int new=0; // 
unsigned int p=1000000000; // 
//
unsigned long long int lnew= 0; //
unsigned long long int lold = (unsigned long long int) old; //
unsigned long long int lp = (unsigned long long int) p; //

printf("unsigned long long int \tunsigned  int \n"); // header 

for ( i=0 ; i<20; i++){
     printf("lnew = %12llu \tnew =  %12u",  lnew, new);
     // check overflow
     // http://stackoverflow.com/questions/2633661/how-to-check-integer-overflow-in-c/
     if ( new < old) printf("    unsigned integer overflow = wrap \n");
                     else printf("\n");          
                       
                        
     // unsigned int
     old=new; // save old value for comparison = overflow check 
     new = old + p ; // simple addition ; new value should be greater then old value 
     // unsigned long long int
     lold=lnew;
     lnew=lold+lp;      
}


return 0; 
}

Wynik :

unsigned long long int 	unsigned  int 
lnew =            0 	new =             0
lnew =   1000000000 	new =    1000000000
lnew =   2000000000 	new =    2000000000
lnew =   3000000000 	new =    3000000000
lnew =   4000000000 	new =    4000000000
lnew =   5000000000 	new =     705032704    unsigned integer overflow = wrap 
lnew =   6000000000 	new =    1705032704
lnew =   7000000000 	new =    2705032704
lnew =   8000000000 	new =    3705032704
lnew =   9000000000 	new =     410065408    unsigned integer overflow = wrap 
lnew =  10000000000 	new =    1410065408
lnew =  11000000000 	new =    2410065408
lnew =  12000000000 	new =    3410065408
lnew =  13000000000 	new =     115098112    unsigned integer overflow = wrap 
lnew =  14000000000 	new =    1115098112
lnew =  15000000000 	new =    2115098112
lnew =  16000000000 	new =    3115098112
lnew =  17000000000 	new =    4115098112
lnew =  18000000000 	new =     820130816    unsigned integer overflow = wrap 
lnew =  19000000000 	new =    1820130816




Zapobieganie
edytuj
  • sprawdzanie danych :[30]
    • przed wykonaniem działań [31][32]
    • po wykonaniu działań ( może być niebezpieczne dla liczb ze znakiem ponieważ niezdefiniowane zachowanie zagraża bezpieczeństwu komputera )
  • zwiększenie limitów poprzez :
    • zmianę typu ( int , long int, long long int )
    • użycie biblioteki o dowolnej precyzji ( np. GMP )


Zapobieganie: wykrywanie możliwego przepełnienia przed wykonaniem działania. Porównaj:

  • kod z scaler
  • kod z c-FAQ[34]

rozmiar

edytuj
/*
Here is a small C program 
that will print out the size in bytes 
of some basic C types on your machine. 


Paul Gribble | Summer 2012
This work is licensed under a Creative Commons Attribution 4.0 International License
http://gribblelab.org/CBootcamp/3_Basic_Types_Operators_And_Expressions.html

gcc b.c -Wall
./a.out
*/


#include <stdio.h>

int main(int argc, char *argv[]) {
        printf("a char is %ld bytes\n", sizeof(char));
        printf("an int is %ld bytes\n", sizeof(int));
        printf("an float is %ld bytes\n", sizeof(float));
        printf("a double is %ld bytes\n", sizeof(double));
        printf("a short int is %ld bytes\n", sizeof(short int));
        printf("a long int is %ld bytes\n", sizeof(long int));
        printf("a long double is %ld bytes\n", sizeof(long double));
        return 0;
}



a char is 1 bytes
an int is 4 bytes
an float is 4 bytes
a double is 8 bytes
a short int is 2 bytes
a long int is 8 bytes
a long double is 16 bytes


Liczba cyfr

edytuj

Liczba cyfr w liczbie zmiennoprzecinkowej [35]

// http://ubuntuforums.org/showthread.php?t=986212
// http://www.cplusplus.com/reference/cfloat/
// gcc d.c -lm -Wall
// ./a.out

#include <stdio.h>
#include <float.h> 

int main(void)
{
    printf("Float can ensure %d decimal places\n", FLT_DIG);
    printf("Double can ensure %d decimal places\n", DBL_DIG);
    printf("Long double can ensure %d decimal places\n", LDBL_DIG);

    return 0;
}

Wynik :

Float can ensure 6 decimal places
Double can ensure 15 decimal places
Long double can ensure 18 decimal places

Liczby subnormalne

edytuj
przybliżenia DBL_MIN i liczby subnormalnej
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) 


eliminacja liczb subnormalnych
edytuj

Ten program generuje liczby subnormale:

/*
https://blogs.oracle.com/d/subnormal-numbers
gcc -O0 f.c 
*/
#include <stdio.h>
void main()
{
  double d=1.0;
  while (d>0) {printf("%e\\n",d); d=d/2.0;}
}


wynik:

1.000000e+00
5.000000e-01
2.500000e-01
1.250000e-01
6.250000e-02
3.125000e-02
...
3.162020e-322
1.581010e-322
7.905050e-323
3.952525e-323
1.976263e-323
9.881313e-324
4.940656e-324

Jeśli jednak skompilujemy go z opcję:

  gcc -O0 -ffast-math f.c

to otrzymamy:

...
3.560118e-307
1.780059e-307
8.900295e-308
4.450148e-308
2.225074e-308

Liczby subnormalne są zaokrąglane do zera.

Jaki wpływ na obliczenia mają liczby subnormalne?
edytuj
  • wydłużają czas obliczeń[36]

część ułamkowa

edytuj

Za pomocą:[37]

  • funkcji modf
  • konwersji int
double frac = r - (int)r ;

Błędy w obliczeniach numerycznych

edytuj

Typy:[38][39]

  • wg etapu operacji :[40]
    • blędne dane wejściowe : niezgodne z oczekiwanym typem
    • dane wejściowe powodują błąd rezultatu
  • wg rodzaju operacji


Efekt:

  • zachowanie niezdefiniowane : ( ang. undefined behavior = UB)[51] [52]

Dzielenie przez zero

edytuj

Dzielenie przez zero[53]

Kiedy dzielnik ma wartość zero w wyrażeniu (zwykle przez pomyłkę) to następuje awaria programu (nieprawidłowe zakończenie an. crash). Nieprawidłowe zakończenie może być poważnym problemem w przypadku oprogramowania o krytycznym znaczeniu dla życia.

Zapobiegać temu można przez:

  • kontrole if-else
  • obsługę wyjątków


Mnożenie

edytuj
#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;
}

Kompilacja i uruchomienie:

  gcc t.c -Wall
  ./a.out

Wynik:



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

Porównywanie

edytuj

Sprawdźmy czy liczba x jest równa zero :

if (x==0.0)


Czy takie porównanie jest bezpieczne dla liczb zmiennoprzecinkowych ?[54]


//  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;
}

Wynik :


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


Jak powinno się porównywać liczby zmienno przecinkowe ?[55]

  • wartość bezwględna róznicy : if( abs(a-b) < epsilon) // wrong - don't do this[56]
  • if( abs((a-b)/b) < epsilon ) // still not right!
  • wartości graniczne
  • stałe [57]
 if (fpclassify(x) == FP_ZERO )

lub

 if (x == FP_ZERO)

Wartości służące do testo wania porównań :

  • wg Michael Borgwardt[58]

Sumowanie

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

Biblioteki matematyczne

edytuj
 
Dependency Graph for tgmath.h

float.h

edytuj

opis

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[60], czyli na końcu linii :[61]

Flaga lm jest zależna od środowiska. Na przykład, w systemie Windows, nie jest to wymagana, ale jest to wymagana w systemach opartych na systemie UNIX.

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.


Wersje

Stałe matematyczne: pi, e ...

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


Typy

  • _Complex
  • complex.h
  • biblioteki:


różnica pomiędzy _complex a complex

edytuj
  • _Complex jest słowem kluczowym, możemy go używać bez dyrektywy include dołączającej plik nagłówkowy complex.h
  • complex jest makrem z pliku nagłówkowego complex.h
float _Complex     
double _Complex                  
long double _Complex


 #include <complex.h>
 float complex
double complex
long double complex

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

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>

która powoduje dołączenie standardowej biblioteki obsługującej liczny zespolenie

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

Więcej:

// https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
// program by user870774
#include <stdio.h>      /* Standard Library of Input and Output */
#include <complex.h>    /* Standard Library of Complex Numbers */

int main() {

    double complex z1 = 1.0 + 3.0 * I;
    double complex z2 = 1.0 - 4.0 * I;

    printf("Working with complex numbers:\n\v");

    printf("Starting values: Z1 = %.2f + %.2fi\tZ2 = %.2f %+.2fi\n", creal(z1), cimag(z1), creal(z2), cimag(z2));

    double complex sum = z1 + z2;
    printf("The sum: Z1 + Z2 = %.2f %+.2fi\n", creal(sum), cimag(sum));

    double complex difference = z1 - z2;
    printf("The difference: Z1 - Z2 = %.2f %+.2fi\n", creal(difference), cimag(difference));

    double complex product = z1 * z2;
    printf("The product: Z1 x Z2 = %.2f %+.2fi\n", creal(product), cimag(product));

    double complex quotient = z1 / z2;
    printf("The quotient: Z1 / Z2 = %.2f %+.2fi\n", creal(quotient), cimag(quotient));

    double complex conjugate = conj(z1);
    printf("The conjugate of Z1 = %.2f %+.2fi\n", creal(conjugate), cimag(conjugate));

    return 0;
}

Opis MPC

Dodatkowe

edytuj

Zobacz również

edytuj

Źródła

edytuj
  1. The Rational Number Class in C by Peter Burden
  2. Rational Arithmetic by R. Sedgewick
  3. Where is the itoa function in Linux?
  4. gcc - Binary-constants
  5. Code Replay : C represent int in base 2
  6. what-every-computer-programmer-should by Josh Haberman
  7. Metody numeryczne - autorzy : Piotr Krzyżanowski i Leszek Plaskota — Uniwersytet Warszawski, Wydział Matematyki, Informatyki i Mechaniki
  8. How to get memory representation-double
  9. IEEE-754 Floating Point Converter
  10. dumpfp: A Tool to Inspect Floating-Point Numbers by Joshua Haberman
  11. subethasoftware : converting-two-8-bit-values-to-one-16-bit-value-in-c
  12. subethasoftware: splitting-a-16-bit-value-to-two-8-bit-values-in-c
  13. Math to Code by cameron smith
  14. developing-mathematical-software-in-c by Fredrik Johansson
  15. printing-algebraic-numbers by Fredrik Johansson
  16. numeryczne - Wydziału MIM UW
  17. Cpp Core Guidelines : arithmetic
  18. i ból obliczeń numerycznych -Piotr Krzyżanowski
  19. Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg
  20. Two disasters caused by computer arithmetic errors
  21. Praktyczne wyznaczanie precyzji arytmetyki - autorzy : Piotr Krzyżanowski i Leszek Plaskota — Uniwersytet Warszawski, Wydział Matematyki, Informatyki i Mechaniki
  22. stackoverflow question : why-does-gcc-report-a-floating-point-exception-when-i-execute-1-0
  23. stackoverflow question: catch-floating-point-exceptions-using-a-compiler-option-with-c
  24. Point Representation - Basics from : geeksforgeeks.org
  25. Tutorial on Data Representation by Chua Hock-Chuan
  26. Przekroczenie zakresu liczb całkowitych w wikipedii
  27. [http://www.fefe.de/intof.html%7CCatching Integer Overflows in C ]
  28. Guide to Undefined Behavior in C and C++, Part 1 by John Regehr
  29. stackoverflow question: function-abs-returning-negative-number-in-c
  30. Testing for Integer Overflow in C and C++ by Josh Haberman
  31. comp.lang.c FAQ list · Question 20.6b : How can I ensure that integer arithmetic doesn't overflow ?
  32. GCC : Built-in Functions to Perform Arithmetic with Overflow Checking
  33. gnulib : Checking-Integer-Overflow
  34. c-faq: three functions for ``careful addition, subtraction, and multiplication
  35. Stackoverflow : Counting digits in a float
  36. O N SUBNORMAL FLOATING POINT AND ABNORMAL TIMING by Marc Andrysco, David Kohlbrenner, Keaton Mowery, Ranjit Jhala, Sorin Lerner, and Hovav Shacham
  37. stackoverflow question : extract-decimal-part-from-a-floating-point-number-in-c
  38. INTRODUCTION TO NUMERICAL ANALYSIS WITH C PROGRAMS by Attila Mate
  39. Basic Issues in Floating Point Arithmetic and Error Analysis by Jim Demmel
  40. securecoding.cert.org : Prevent+or+detect+domain+and+range+errors+in+math+functions
  41. Floating point inaccuracy examples
  42. Catastrophic Cancellation: The Pitfalls of Floating Point Arithmetic by Graham Markall
  43. math.stackexchange question: is-this-characteristic-of-tent-map-usually-observed
  44. BŁĘDY PRZETWARZANIA NUMERYCZNEGO, Maciej Patan, UW Zielonogórski
  45. dangerous is it to compare floating point values?
  46. it possible to get 0 by subtracting two unequal floating point numbers?
  47. How to solve quadratic equations numerically by FLORIAN DANG
  48. Double Rounding Errors in Floating-Point Conversions By Rick Regan (Published August 11th, 2010)
  49. stackoverflow question when-does-underflow-occur
  50. stackoverflow question: how-to-define-underflow-for-an-implementationieee754-which-support-subnormal-n
  51. Undefined_behavior w ang. wikipedii
  52. undefined-behavior-in-c-and-cplusplus-programs by Nayuki
  53. delftstack :infinity-in-c by Abdul Mateen Oct-26, 2022
  54. stackoverflow question: c-floating-point-zero-comparison
  55. Comparing Floating-Point Numbers Is Tricky by Matt Kline
  56. | floating-point-gui.de : comparison/
  57. stackoverflow question : float-double-equality-with-exact-zero
  58. Michael Borgwardt : Nearly Equals Test in java
  59. numerically-stable-law-of-cosines by Nayuki
  60. Stackoverflow : Undefined reference to `pow' and `floor'
  61. I'm using math.h and the library link -lm, but “undefined reference to `pow'” still happening
  62. arb library