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\left(0\right)=\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 :

#!/usr/bin/env python3
# file: dct.py
# (C) 2019 by marnout à free pt fr  Released under the GPL
 
import mpmath as mp
 
class Dct:
   
   def __init__(self, data='data'): 
      mp.mp.prec = 80
      self.u = list()
      with open(data) as f:
         for data in f:
            self.u.append(mp.mpf(data))
      self.N = len(self.u)
 
   @mp.workprec(100)
   def _dct(self):
      self.dct = list()
      # dct type II avec normalisation
      # dct(k) = 2γ(k)}/N*Σ_n=0^n=N-1 u[n].cos((2n+1)kπ/2N)
      # γ(0) = 1/sqrt(2) et 1 autrement
      for k in range(self.N):
         s = mp.mpf(0) # sum
         for n in range(self.N):
            t = mp.fmul(mp.mpf((2*n + 1)*k), mp.pi()) # t = (2n + 1)kπ
            t = mp.fdiv(t, self.N*2) # t = (2n + 1)kπ/2N
            t = mp.cos(t) # t = cos((2n + 1)kπ/2N)
            t = mp.fmul(self.u[n], t) # t = u[n]cos((2n + 1)kπ/2N)
            s = mp.fadd(s, t) # Σ
         s = mp.fmul(s, 2)
         s = mp.fdiv(s, self.N)
         if k == 0:
            s = mp.fdiv(s, mp.sqrt(2)) # orthogonalisation
         self.dct.append(s)
 
   def idct(self, n):
      s = mp.mpf(0) # sum
      # idct(n) = Σ_k=0^k=N-1 γ(k)dct(k)cos((2n+1)kπ/2N))
      for k in range(self.N):
         t = mp.fmul(mp.mpf((2*n + 1)*k), mp.pi()) # t = (2n + 1)kπ
         t = mp.fdiv(t, 2*self.N) # t = (2n + 1)kπ/2N
         t = mp.cos(t) # t = cos((2n + 1)kπ/2N)
         t = mp.fmul(self.dct[k], t) # t = dct(k)cos((2n+1)kπ/2N)
         if k == 0: t =  mp.fdiv(t, mp.sqrt(2)) # Σ
         s = mp.fadd(s, t)
         # print("k:", k, " n:", n, ' s:', s, file=f)
      return s
 
 
if __name__ == '__main__':
   
   print(f"\n{__file__} - Mourad Arnout - LGPL\n")
 
   d = Dct()
   print("Precision :", mp.mp.prec, '\n')
   
   print("dct")
   d._dct()
   for k in range(d.N): 
      print("dct(%2d) = %24.6f" % (k, d.dct[k]))
 
   print("idct")  
   for k in range(d.N): 
      print("u[%2d] = %16.0f" % (k,  d.u[k]), end=' ; ')
      print("idct(%2d) = %24.6f" % (k, d.idct(k)))
      
 
   print(f"\nfin {__file__}")

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

/* 
dct.c
Released under the GPL (C) 2019 marnout à free pt fr
*/
#include <stdio.h>
#include <string.h>
#include <mpfr.h>
 
#define RND MPFR_RNDN
 
const mpfr_prec_t prec = 80; // précision des nombres
const size_t N = 16; // nombre de données 
const size_t wdsize = 20; // taille d'un mot
 
// int main(int argc, char *argv[]) 
int main(void) 
{
   mpfr_set_default_prec(prec); // précision des nombres
   mpfr_t u[N], dct[N], s, t; 
   mpfr_inits(s, t, (mpfr_ptr) 0);
   // read data
   char data[wdsize];
   FILE *fp = fopen("data", "r");
   for(size_t k=0; k<N; k++) {
      mpfr_init(u[k]);
      if(fgets(data, wdsize, fp) == NULL) break;
      mpfr_set_str(u[k], data, 0, RND);
   }
   fclose(fp);
   // dct(k) = 2γ(k)}/N*Σ_n=0^n=N-1 u[n].cos((2n+1)kπ/2N)
   // γ(0) = 1/sqrt(2) et 1 autrement
   mpfr_t sqrt2; mpfr_init(sqrt2); mpfr_sqrt_ui(sqrt2, 2, RND);
   printf("DCT\n");
   for(size_t k=0; k<N; k++) {
      mpfr_init(dct[k]);
      mpfr_set_zero(s, 0);
      for(size_t n=0; n<N; n++) { // calcul de s
         mpfr_const_pi(t, RND); // t = π
         mpfr_mul_ui(t, t, (2*n+1)*k, RND); // t = (2n+1)kπ
         mpfr_div_ui(t, t, 2*N, RND); // t = (2n+1)kπ/2N
         mpfr_cos(t, t, RND); // t = cos((2n+1)kπ/2N)
         mpfr_mul(t, u[n], t, RND); // t = u[n]cos((2n+1)kπ/2N)
         mpfr_add(s, s, t, RND); // Σ
      }
      mpfr_mul_ui(s, s, 2, RND);
      mpfr_div_ui(s, s, N, RND);
      if(k == 0) mpfr_div(s, s, sqrt2, RND);
      mpfr_set(dct[k], s, RND);
      mpfr_printf("DCT(%2u) = %25Rf\n", k, dct[k]);
   } // dct
   printf("\nIDCT\n");
   // idct(n) = Σ_k=0^k=N-1 γ(k)dct(k)cos((2n+1)kπ/2N))
   for(size_t n=0; n<N; n++) {
      mpfr_set_zero(s, 0); // s = 0
      for(size_t k=0; k<N; k++) { // calcul de s
         mpfr_const_pi(t, RND); // t = π
         mpfr_mul_ui(t, t, k*(2*n+1), RND); // t = (2n+1)kπ
         mpfr_div_ui(t, t, 2*N, RND); // t = (2n+1)kπ/2N
         mpfr_cos(t, t, RND); // t = cos((2n+1)kπ/2N)
         mpfr_mul(t, t, dct[k], RND); // t = dct(k)cos((2n+1)kπ/2N)
         if(k == 0) mpfr_div(t, t, sqrt2, RND);
         mpfr_add(s, s, t, RND);
      }
      mpfr_round(s, s);
      mpfr_printf("u[%2d] = %16.0Rf", n, u[n]);
      mpfr_printf(" ; IDCT(%2d) = %16.0Rf\n", n, s);
      
   } // idct
   mpfr_clears(s, u, t, (mpfr_ptr) 0);
   for(size_t k=0; k<N; k++) mpfr_clear(dct[k]);
   mpfr_free_cache();
   return 0;
}

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