jeudi 28 mars 2019

Calculs en multiple précision

Pense-bête pour les calculs en multiple précision en python (module mpmath) et en C (GNU MPFR).

Installation et premiers pas

Python

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.

Langage C

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

Précision et arrondi

Python

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.

Langage C

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

Exemple : DCT (discrete cosine transform)

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)$.

Python

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)$.

Langage C

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)$


Réalisé avec Qlam - LGPL