C/Zaawansowane operacje matematyczne
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
edytujLiczby 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- całkowite
- wymierne
- z pomocą typu mpq_t z biblioteki GMP ( GNU Multiple Precision Arithmetic Library )
- z użyciem struktur [1][2]
- zmiennoprzecinkowe
- zespolone
- kwaterniony
Dodatkowo rodzaj liczby definiują :
- specyfikatory
- podstawa systemu liczbowego
- binarne
- ósemkowe (oktalne)
- dziesiętne
- szesnastkowe
Jak widać nie ma :
- Liczb niewymiernych ( chyba że korzystamy z obliczeń symbolicznych, np. program Maxima CAS )
liczby binarne
edytujSposoby :
- itoa
- snprintf [3]
- Binarna stała całkowita z użyciem rozszerzenia gcc
- za pomoc operacji na bitach, np. shift
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)
edytujLiczba 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- How Computers Use Numbers
- Floating Point Math by Erik Wiffin
- float exposed by Bartosz Ciechanowski
- The Theory Behind How Computers Subtract by Rutvik
Liczby rzeczywisty
edytujW 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ż:
- dumpfp[10]
Liczby całkowite
edytujLiczba 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
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:
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
edytujsumowanie
edytujW 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)
edytujW 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
edytujObliczenia numeryczne
edytujObliczenia numeryczne[16] [17]są to obliczenia na liczbach. Ich przeciwieństwem są obliczenia symboliczne wykonywane na symbolach ( zobacz Maxima CAS ) [18]
Uwaga!
|
Epsilon maszynowy
edytujEpsilon 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ń.
Uwaga!
|
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
edytujWyjatki
- 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ń
edytujzmiennoprzecinkowych
edytujDefinicje 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)
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
edytujPrzekroczenie 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]
- 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
edytujLiczba 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
edytujprzybliżenia DBL_MIN i liczby subnormalnej
edytujKorzystają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
edytujTen 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
edytujZa pomocą:[37]
- funkcji modf
- konwersji int
double frac = r - (int)r ;
Błędy w obliczeniach numerycznych
edytuj- wg etapu operacji :[40]
- blędne dane wejściowe : niezgodne z oczekiwanym typem
- dane wejściowe powodują błąd rezultatu
- wg rodzaju operacji
Efekt:
Dzielenie przez zero
edytujDzielenie 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
edytujSprawdź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
edytujNa 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
edytujUtrata 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
edytujfloat.h
edytujmath.h
edytujAby 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
- standardowa: math.h
- openlibm: a high quality, portable, standalone C mathematical library (libm)
- crlibm: a library of correct rounding elementary functions in double precision ( correctly rounded libm)
- tgmath.h
Stałe matematyczne: pi, e ...
edytujW 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:
- mpc
- arb [62]
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
edytujOperacje na liczbach zespolonych są częścią uaktualnionego standardu języka C o nazwie C99, który jest obsługiwany jedynie przez część kompilatorów
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>
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;
}
MPC
edytujDodatkowe
edytujZobacz również
edytujŹródła
edytuj- ↑ The Rational Number Class in C by Peter Burden
- ↑ Rational Arithmetic by R. Sedgewick
- ↑ Where is the itoa function in Linux?
- ↑ gcc - Binary-constants
- ↑ Code Replay : C represent int in base 2
- ↑ what-every-computer-programmer-should by Josh Haberman
- ↑ Metody numeryczne - autorzy : Piotr Krzyżanowski i Leszek Plaskota — Uniwersytet Warszawski, Wydział Matematyki, Informatyki i Mechaniki
- ↑ How to get memory representation-double
- ↑ IEEE-754 Floating Point Converter
- ↑ dumpfp: A Tool to Inspect Floating-Point Numbers by Joshua Haberman
- ↑ subethasoftware : converting-two-8-bit-values-to-one-16-bit-value-in-c
- ↑ subethasoftware: splitting-a-16-bit-value-to-two-8-bit-values-in-c
- ↑ Math to Code by cameron smith
- ↑ developing-mathematical-software-in-c by Fredrik Johansson
- ↑ printing-algebraic-numbers by Fredrik Johansson
- ↑ numeryczne - Wydziału MIM UW
- ↑ Cpp Core Guidelines : arithmetic
- ↑ i ból obliczeń numerycznych -Piotr Krzyżanowski
- ↑ Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg
- ↑ Two disasters caused by computer arithmetic errors
- ↑ Praktyczne wyznaczanie precyzji arytmetyki - autorzy : Piotr Krzyżanowski i Leszek Plaskota — Uniwersytet Warszawski, Wydział Matematyki, Informatyki i Mechaniki
- ↑ stackoverflow question : why-does-gcc-report-a-floating-point-exception-when-i-execute-1-0
- ↑ stackoverflow question: catch-floating-point-exceptions-using-a-compiler-option-with-c
- ↑ Point Representation - Basics from : geeksforgeeks.org
- ↑ Tutorial on Data Representation by Chua Hock-Chuan
- ↑ Przekroczenie zakresu liczb całkowitych w wikipedii
- ↑ [http://www.fefe.de/intof.html%7CCatching Integer Overflows in C ]
- ↑ Guide to Undefined Behavior in C and C++, Part 1 by John Regehr
- ↑ stackoverflow question: function-abs-returning-negative-number-in-c
- ↑ Testing for Integer Overflow in C and C++ by Josh Haberman
- ↑ comp.lang.c FAQ list · Question 20.6b : How can I ensure that integer arithmetic doesn't overflow ?
- ↑ GCC : Built-in Functions to Perform Arithmetic with Overflow Checking
- ↑ gnulib : Checking-Integer-Overflow
- ↑ c-faq: three functions for ``careful addition, subtraction, and multiplication
- ↑ Stackoverflow : Counting digits in a float
- ↑ O N SUBNORMAL FLOATING POINT AND ABNORMAL TIMING by Marc Andrysco, David Kohlbrenner, Keaton Mowery, Ranjit Jhala, Sorin Lerner, and Hovav Shacham
- ↑ stackoverflow question : extract-decimal-part-from-a-floating-point-number-in-c
- ↑ INTRODUCTION TO NUMERICAL ANALYSIS WITH C PROGRAMS by Attila Mate
- ↑ Basic Issues in Floating Point Arithmetic and Error Analysis by Jim Demmel
- ↑ securecoding.cert.org : Prevent+or+detect+domain+and+range+errors+in+math+functions
- ↑ Floating point inaccuracy examples
- ↑ Catastrophic Cancellation: The Pitfalls of Floating Point Arithmetic by Graham Markall
- ↑ math.stackexchange question: is-this-characteristic-of-tent-map-usually-observed
- ↑ BŁĘDY PRZETWARZANIA NUMERYCZNEGO, Maciej Patan, UW Zielonogórski
- ↑ dangerous is it to compare floating point values?
- ↑ it possible to get 0 by subtracting two unequal floating point numbers?
- ↑ How to solve quadratic equations numerically by FLORIAN DANG
- ↑ Double Rounding Errors in Floating-Point Conversions By Rick Regan (Published August 11th, 2010)
- ↑ stackoverflow question when-does-underflow-occur
- ↑ stackoverflow question: how-to-define-underflow-for-an-implementationieee754-which-support-subnormal-n
- ↑ Undefined_behavior w ang. wikipedii
- ↑ undefined-behavior-in-c-and-cplusplus-programs by Nayuki
- ↑ delftstack :infinity-in-c by Abdul Mateen Oct-26, 2022
- ↑ stackoverflow question: c-floating-point-zero-comparison
- ↑ Comparing Floating-Point Numbers Is Tricky by Matt Kline
- ↑ | floating-point-gui.de : comparison/
- ↑ stackoverflow question : float-double-equality-with-exact-zero
- ↑ Michael Borgwardt : Nearly Equals Test in java
- ↑ numerically-stable-law-of-cosines by Nayuki
- ↑ Stackoverflow : Undefined reference to `pow' and `floor'
- ↑ I'm using math.h and the library link -lm, but “undefined reference to `pow'” still happening
- ↑ arb library