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

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

En 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

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 et mpmath

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):
    26	            t = mp.fmul(mp.mpf((2*n + 1)*k), mp.pi()) # t = (2n + 1)kπ
    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)
    32	         s = mp.fdiv(s, self.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	         # print("k:", k, " n:", n, ' s:', s, file=f)
    48	      return s
    49	
    50	
    51	if __name__ == '__main__':
    52	   
    53	   print(f"\n{__file__} - Mourad Arnout - LGPL\n")
    54	
    55	   d = Dct()
    56	   print("Precision :", mp.mp.prec, '\n')
    57	   
    58	   print("dct")
    59	   d._dct()
    60	   for k in range(d.N): 
    61	      print("dct(%2d) = %24.6f" % (k, d.dct[k]))
    62	
    63	   print("idct")  
    64	   for k in range(d.N): 
    65	      print("u[%2d] = %16.0f" % (k,  d.u[k]), end=' ; ')
    66	      print("idct(%2d) = %24.6f" % (k, d.idct(k)))
    67	      
    68	
    69	   print(f"\nfin {__file__}")
    70	

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 et MPFR

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