jeudi 28 mars 2019
Pense-bête pour les calculs en multiple précision en python (module mpmath) et en C (GNU MPFR).
En python il faut installer le paquet python-mpmath.
Utilisation :
>>> from mpmath import *
>>> mpf(2)**mpf(0.5) # racine carrée de 2
mpf('1.4142135623730951')
>>> mp.prec = 128 # precision de 128 bits, soit 38 chiffres en base dix
>>> mpf(2)**mpf(0.5)
mpf('1.414213562373095048801688724209698078569')
mpf : muliple précision float
mp.prec : précision en bits
La documentation se trouve à mpmath.org.
En C on installe le paquet mpfr qui est maintenu par l'Inria.
Utilisation \n
#include <stdio.h>
#include <mpfr.h>
//
void main(void)
{
mpfr_t x, a; // declare x et a
mpfr_init2(x, 128); // initialise x avec 128 bits
mpfr_init2(a, 128); // initialise a
mpfr_set_d(a, -1.0, MPFR_RNDN); // affecte -1 à a (a := -1)
mpfr_acos(x, a, MPFR_RNDN); // x = acos(a)
printf("x: "); // affiche x
pfr_out_str(stdout, 10, 0, x, MPFR_RNDN);
mpfr_const_pi(x, MPFR_RNDN); // utilise la constante pi de mpfr
mpfr_printf("\nx: %.45Rf", x); // analogue à printf de C
mpfr_clear(x); // libere la mémoire
mpfr_clear(a);
}
La compilation se fait avec
gcc test.c -o test -lmpfr -lgmp
On obtient la sortie
x: 3.141592653589793238462643383279502884195 x: 3.141592653589793238462643383279502884195286358
mpfr_t : multiple precision floating-point reliably type
MPFR_RNDN : mode d'arrondi, ici le plus proche (nearest)
mpfr_out_str : sortie sur le flux (output on stream)
mpfr_printf est compatible avec printf de C. Elle ajoute entre autres le type R pour afficher les mpfr_t et le type P pour les mpfr_prec_t.
Documentation : GNU MPFR
On a déjà vu mp.prec, la précision en bits (binaire)
Il y a aussi mp.dps (decimal places) pour la précision décimale. Le rapport entre les deux est de 3,33.
On peut modifier temporairement la précision qui est par défaut mp.prec=53 avec mp.prec += 75 puis mp.prec -= 75.
Il y a aussi workprec sous forme de directive ou de décorateur. Voir Precision management dans la doc officielle.
Nous avons rencontré mpfr_init2 pour l'initialisation d'une variable. On peut utiliser mpfr_init(mpfr_t x) qui utilise la précision par défaut (53 bits) ou celle fixée par le programmeur à l'aide de void mpfr_set_default_prec(mpfr_prec_t prec). Cf. Iniitialisation Functions.
De plus mpfr utilise plusieurs modes d'arrondis. Voir Rounding Modes
C'est un succédané de la transformation de Fourier. Cette variante n'utilise que les réels et la fonction cosinus. Son intérêt est qu'elle est mieux adaptée à des fonctions non périodiques.
Soit $f$ une fonction définie sur $[0, N-1]$,
pour $k\in \{0, 1, ... , N-1\}$, on note $u_k=f(k)$ .
La tranformée en cosinus de f est
$$dct\left(k\right)=~\frac{2c(k)}{N}~\sum_{n=0}^{N-1}{u_n .cos\frac{(2n+1)k\pi}{2N}}$$
où $c(0)=\dfrac{1}{\sqrt 2}$ et $c(k)=1$ pour $k\ne 0$
La transformation inverse sera obtenue par
$$idct\left(n\right)=\sum_{k=0}^{N-1}{c(k).dct(k).cos\frac{(2n+1)k\pi}{2N}}$$
Remarque : d'autres variantes existent voir Transformée en cosinus discrète sur Wikipedia ou Interpolation et approximation etc.
Remarque : Cette transformation est utilisée par la compression JPEG.
Pour illustrer notre exemple on prendra la suite de nombre enregistrée dan le fichier data :
2115088 17592193908736 35188802846720 2199092987904 423346336448512 68720553984 13194677059584 4362600484 277293827328 275146342816 146163113986 6597080252672 51842646016 17635269935360 70369012875784 140738562164800
Nous voudrions afficher les $dct(k)$ puis les $idct(n)$.
Script python avec le module mpmath :
1 #!/usr/bin/env python3 2 # file: dct.py 3 # (C) 2019 by marnout à free pt fr Released under the GPL 4 5 import mpmath as mp 6 7 class Dct: 8 9 def __init__(self, data='data'): 10 mp.mp.prec = 80 11 self.u = list() 12 with open(data) as f: 13 for data in f: 14 self.u.append(mp.mpf(data)) 15 self.N = len(self.u) 16 17 @mp.workprec(100) 18 def _dct(self): 19 self.dct = list() 20 # dct type II avec normalisation 21 # dct(k) = 2γ(k)}/N*Σ_n=0^n=N-1 u[n].cos((2n+1)kπ/2N) 22 # γ(0) = 1/sqrt(2) et 1 autrement 23 for k in range(self.N): 24 s = mp.mpf(0) # sum 25 for n in range(self.N): # t = (2n + 1)kπ 26 t = mp.fmul(mp.mpf((2*n + 1)*k), mp.pi()) 27 t = mp.fdiv(t, self.N*2) # t = (2n + 1)kπ/2N 28 t = mp.cos(t) # t = cos((2n+1)kπ/2N) 29 t = mp.fmul(self.u[n], t) # t = u[n]cos((2n+1)kπ/2N) 30 s = mp.fadd(s, t) # Σ 31 s = mp.fmul(s, 2) # s = 2Σ 32 s = mp.fdiv(s, self.N) # s = 2/NΣ 33 if k == 0: 34 s = mp.fdiv(s, mp.sqrt(2)) # orthogonalisation 35 self.dct.append(s) 36 37 def idct(self, n): 38 s = mp.mpf(0) # sum 39 # idct(n) = Σ_k=0^k=N-1 γ(k)dct(k)cos((2n+1)kπ/2N)) 40 for k in range(self.N): 41 t = mp.fmul(mp.mpf((2*n + 1)*k), mp.pi()) # t = (2n+1)kπ 42 t = mp.fdiv(t, 2*self.N) # t = (2n+1)kπ/2N 43 t = mp.cos(t) # t = cos((2n+1)kπ/2N) 44 t = mp.fmul(self.dct[k], t) # t = dct(k)cos((2n+1)kπ/2N) 45 if k == 0: t = mp.fdiv(t, mp.sqrt(2)) # Σ 46 s = mp.fadd(s, t) 47 return s 48 49 50 if __name__ == '__main__': 51 52 print(f"∖n{__file__} - MArnout - LGPL∖n") 53 54 d = Dct() 55 print("Precision :", mp.mp.prec, '∖n') 56 57 print("dct") 58 d._dct() 59 for k in range(d.N): 60 print("dct(\%2d) = \%24.6f" \% (k, d.dct[k])) 61 62 print("idct") 63 for k in range(d.N): 64 print("u[\%2d] = \%16.0f" \% (k, d.u[k]), end=' ; ') 65 print("idct(\%2d) = \%24.6f" \% (k, d.idct(k))) 66 67 68 print(f"∖nfin {__file__}") 69
Voici dct.py en téléchargement.
Résultat du script :
Precision : 80 dct dct( 0) = 64318835839289.390625 dct( 1) = 11830695203688.791016 dct( 2) = 18185603126586.210938 dct( 3) = -68102748744173.234375 dct( 4) = -31311371966301.996094 dct( 5) = -31651701025609.765625 dct( 6) = 36033499476967.250000 dct( 7) = 37874797111808.484375 dct( 8) = 37012142079090.171875 dct( 9) = -10074053223196.046875 dct(10) = -42547016635737.375000 dct(11) = -49927356783822.343750 dct(12) = -19516995321046.324219 dct(13) = 25438436714700.421875 dct(14) = 56154252525608.046875 dct(15) = 42911581796127.429688 idct u[ 0] = 2115088 ; idct( 0) = 2115088.000000 u[ 1] = 17592193908736 ; idct( 1) = 17592193908736.000000 u[ 2] = 35188802846720 ; idct( 2) = 35188802846720.000000 u[ 3] = 2199092987904 ; idct( 3) = 2199092987904.000000 u[ 4] = 423346336448512 ; idct( 4) = 423346336448512.000000 u[ 5] = 68720553984 ; idct( 5) = 68720553984.000000 u[ 6] = 13194677059584 ; idct( 6) = 13194677059584.000000 u[ 7] = 4362600484 ; idct( 7) = 4362600484.000000 u[ 8] = 277293827328 ; idct( 8) = 277293827328.000000 u[ 9] = 275146342816 ; idct( 9) = 275146342816.000000 u[10] = 146163113986 ; idct(10) = 146163113986.000000 u[11] = 6597080252672 ; idct(11) = 6597080252672.000000 u[12] = 51842646016 ; idct(12) = 51842646016.000000 u[13] = 17635269935360 ; idct(13) = 17635269935360.000000 u[14] = 70369012875784 ; idct(14) = 70369012875784.000000 u[15] = 140738562164800 ; idct(15) = 140738562164800.000000
On voit bien que $u_n=idct(n)$.
Le même programme en C utilisant la librairie MPFR
1 /* 2 dct.c 3 Released under the GPL (C) 2019 marnout à free pt fr 4 */ 5 #include <stdio.h> 6 #include <string.h> 7 #include <mpfr.h> 8 9 #define RND MPFR_RNDN 10 11 const mpfr_prec_t prec = 80; // précision des nombres 12 const size_t N = 16; // nombre de données 13 const size_t wdsize = 20; // taille d'un mot 14 15 // int main(int argc, char *argv[]) 16 int main(void) 17 { 18 mpfr_set_default_prec(prec); // précision des nombres 19 mpfr_t u[N], dct[N], s, t; 20 mpfr_inits(s, t, (mpfr_ptr) 0); 21 // read data 22 char data[wdsize]; 23 FILE *fp = fopen("data", "r"); 24 for(size_t k=0; k<N; k++) { 25 mpfr_init(u[k]); 26 if(fgets(data, wdsize, fp) == NULL) break; 27 mpfr_set_str(u[k], data, 0, RND); 28 } 29 fclose(fp); 30 // dct(k) = 2γ(k)}/N*Σ_n=0^n=N-1 u[n].cos((2n+1)kπ/2N) 31 // γ(0) = 1/sqrt(2) et 1 autrement 32 mpfr_t sqrt2; mpfr_init(sqrt2); mpfr_sqrt_ui(sqrt2, 2, RND); 33 printf("DCT∖n"); 34 for(size_t k=0; k<N; k++) { 35 mpfr_init(dct[k]); 36 mpfr_set_zero(s, 0); 37 for(size_t n=0; n<N; n++) { // calcul de s 38 mpfr_const_pi(t, RND); // t = π 39 mpfr_mul_ui(t, t, (2*n+1)*k, RND); // t = (2n+1)kπ 40 mpfr_div_ui(t, t, 2*N, RND); // t = (2n+1)kπ/2N 41 mpfr_cos(t, t, RND); // t = cos((2n+1)kπ/2N) 42 mpfr_mul(t, u[n], t, RND); // t = u[n]cos((2n+1)kπ/2N) 43 mpfr_add(s, s, t, RND); // Σ 44 } 45 mpfr_mul_ui(s, s, 2, RND); 46 mpfr_div_ui(s, s, N, RND); 47 if(k == 0) mpfr_div(s, s, sqrt2, RND); 48 mpfr_set(dct[k], s, RND); 49 mpfr_printf("DCT(\%2u) = \%25Rf∖n", k, dct[k]); 50 } // dct 51 printf("∖nIDCT∖n"); 52 // idct(n) = Σ_k=0^k=N-1 γ(k)dct(k)cos((2n+1)kπ/2N)) 53 for(size_t n=0; n<N; n++) { 54 mpfr_set_zero(s, 0); // s = 0 55 for(size_t k=0; k<N; k++) { // calcul de s 56 mpfr_const_pi(t, RND); // t = π 57 mpfr_mul_ui(t, t, k*(2*n+1), RND); // t = (2n+1)kπ 58 mpfr_div_ui(t, t, 2*N, RND); // t = (2n+1)kπ/2N 59 mpfr_cos(t, t, RND); // t = cos((2n+1)kπ/2N) 60 mpfr_mul(t, t, dct[k], RND); // t = dct(k)cos((2n+1)kπ/2N) 61 if(k == 0) mpfr_div(t, t, sqrt2, RND); 62 mpfr_add(s, s, t, RND); 63 } 64 mpfr_round(s, s); 65 mpfr_printf("u[\%2d] = \%16.0Rf", n, u[n]); 66 mpfr_printf(" ; IDCT(\%2d) = \%16.0Rf∖n", n, s); 67 68 } // idct 69 mpfr_clears(s, u, t, (mpfr_ptr) 0); 70 for(size_t k=0; k<N; k++) mpfr_clear(dct[k]); 71 mpfr_free_cache(); 72 return 0; 73 }
Voici le fichier dct.c en téléchargement.
Le résultat est semblable au script python
$ ./dct DCT DCT( 0) = 64318835839289.392138 DCT( 1) = 11830695203688.791882 DCT( 2) = 18185603126586.209348 DCT( 3) = -68102748744173.236873 DCT( 4) = -31311371966301.994429 DCT( 5) = -31651701025609.764182 DCT( 6) = 36033499476967.251851 DCT( 7) = 37874797111808.486264 DCT( 8) = 37012142079090.175664 DCT( 9) = -10074053223196.046458 DCT(10) = -42547016635737.375221 DCT(11) = -49927356783822.340006 DCT(12) = -19516995321046.322296 DCT(13) = 25438436714700.421803 DCT(14) = 56154252525608.050736 DCT(15) = 42911581796127.432834 IDCT u[ 0] = 2115088 ; IDCT( 0) = 2115088 u[ 1] = 17592193908736 ; IDCT( 1) = 17592193908736 u[ 2] = 35188802846720 ; IDCT( 2) = 35188802846720 u[ 3] = 2199092987904 ; IDCT( 3) = 2199092987904 u[ 4] = 423346336448512 ; IDCT( 4) = 423346336448512 u[ 5] = 68720553984 ; IDCT( 5) = 68720553984 u[ 6] = 13194677059584 ; IDCT( 6) = 13194677059584 u[ 7] = 4362600484 ; IDCT( 7) = 4362600484 u[ 8] = 277293827328 ; IDCT( 8) = 277293827328 u[ 9] = 275146342816 ; IDCT( 9) = 275146342816 u[10] = 146163113986 ; IDCT(10) = 146163113986 u[11] = 6597080252672 ; IDCT(11) = 6597080252672 u[12] = 51842646016 ; IDCT(12) = 51842646016 u[13] = 17635269935360 ; IDCT(13) = 17635269935360 u[14] = 70369012875784 ; IDCT(14) = 70369012875784 u[15] = 140738562164800 ; IDCT(15) = 140738562164800
On vérifie aussi $u_k=IDCT(k)$