Calcul de sommes

Des expériences

Le calcul de sommes est un favori de l'initiation à l'algo au lycée.
Prenons la somme partielle de la série harmonique :

$$
\sum_{k=1}^{k=n}\frac{1}{k}
$$

Rien de plus simple à programmer :

Pseudo-langage
S <- 0
pour k de 1 jusque n faire
   S <- S + 1/k
fpour
retourne S

Par exemple, on fait un petit programme en python :

Python
def somme1(n):
    S = 0
    for k in range(1,n+1):
        S += 1.0 / k
    return S

On met bien 1.0 et non pas 1 pour faire comprendre qu'on travaille avec des nombres à virgule flottante et non pas des entiers sinon $1/k$ pour $k > 1$ va renvoyer 0 : c'est le quotient de la division euclidienne...

Python
In [1]: somme1(100000)
Out[1]: 12.090146129863335

Oui, mais c'est bon ou pas ? Il n'y a pas de raison que cela soit faux : la somme est finie.

Bon, allez, on essaie avec XCAS

Giac-XCAS
somme1(n):={
  local k,S;
  S:=0;
  pour k de 1 jusque n faire
      S := S + 1.0/k ;
  fpour;
  retourne S
}:;

et on lance le même calcul:

Giac-XCAS
12.0901461272

Tiens! Ce n'est pas pareil.

Bon, allez, on voit ce que ça donne sur Haskell : on utilise sum qui somme les éléments d'une liste et map qui apllique une fonction sur les éléments d'une liste.

Haskell
Prelude> sum (map (\ k -> 1 / k) [1..100000])
12.090146129863335

Ah! Comme python


En Ocaml :

OCaml
utop[1]> BatList.fold_left (+.) 0. [? BatList: 1.0 /. (float k) | k <- (1 -- 100_000) ?];;
- : float = 12.0901461299

Ou bien :

OCaml
utop[2]> BatEnum.fold (+.) 0. (BatEnum.map (fun k -> 1. /. (float k)) (1 -- 100_000));;
- : float = 12.0901461299

Bon, qui a raison ? On peut revenir à XCAS mais en utilisant ses capacités de calcul formel : on calcule la valeur exacte de la somme sous forme d'une fraction (monstrueuse...) PUIS on effectue une approximation de cette fraction:

Giac-XCAS
1>> evalf(somme(1/k,k,1,100000),18)
// Temps mis pour l'évaluation: 3.19
0.120901461298634280e2

XCAS n'est pas d'accord avec lui-même ?

Allons voir du côté de Sage : c'est la même fonction qu'en python mais on peut travailler avec des rationnels et ensuite demander une approximation avec la méthode n(taille de la mantisse).

Sage
sage: def somme1(n):
        S = 0
        for k in range(1,n+1):
                S += 1 / k
        return S
sage: somme1(100000).n(64)
12.0901461298634279

Allez, encore une idée saugrenue ! Calculons la somme dans l'autre sens : de 100 000 jusqu'à 1...

Python
def somme2(n):
    S = 0
    for k in range(n,0,-1):
        S += 1.0 / k
    return S

donne

Python
In [22]: somme2(100000)
Out[22]: 12.090146129863408

Idem pour Haskell

Haskell
Prelude> sum (map (\ k -> 1 / k) [100000, 99999 .. 1])
12.090146129863408

La somme n'est pas commutative !

"Boris Vian"
Y a quelqu' chose qui cloche là-d'dans
J'y retourne immédiatement

Comment la machine fait-elle des additions ?

Et bien pas comme nous ! Mais elle calcule « bien » quand même...selon ses critères et ses capacités.

Le calcul en virgule flottante est un sujet d'étude en soi mais trop méconnu à la fois par les mathématicien(ne)s (trop près de la machine et trop loin des mathématiques « habituelles »...) et les informaticien(ne)s (trop près des mathématiques...et également trop près de la machine pour beaucoup !).

Pour un petit résumé on pourra se référer à ce polyqui lui-même contient de nombreuses références.

Allez, pour vous prouver que je ne suis pas regardant sur le langage, nous allons illustrer notre propos avec JavaScript (en fait, nous nous fichons du langage, du moment qu'il suit la norme IEEE 754 concernant l'arithmétique des flottants). Il faut tout de même savoir que JavaScript est utile pour faire du Web. Pour ne pas surchauffer les cerveaux des développeurs, les concepteurs de JavaScript ont décidé qu'il n'y aurait
qu'un seul type numérique : les nombres à virgule flottante de type binary 64 suivant la norme IEEE 754. Il n'y a donc pas d'entiers en JavaScript par exemple !

Javascript
js> Math.pow(2,52)
4503599627370496
js> Math.pow(2,52) + 1
4503599627370497
 
js> Math.pow(2,53)
9007199254740992
js> Math.pow(2,53) + 1
9007199254740992
js> Math.pow(2,53) + 2
9007199254740994
js> Math.pow(2,53) + 3
9007199254740996
js> Math.pow(2,53) + 4
9007199254740996
js> Math.pow(2,53) + 5
9007199254740996
 
js> Math.pow(2,55)
36028797018963970
js> Math.pow(2,55) + 1
36028797018963970
js> Math.pow(2,55) + 2
36028797018963970
js> Math.pow(2,55) + 3
36028797018963970
js> Math.pow(2,55) + 4
36028797018963970
js> Math.pow(2,55) + 5
36028797018963976
js> Math.pow(2,55) + 10
36028797018963976
js> Math.pow(2,55) + 12
36028797018963980

Il faut avoir en tête que la machine ne travaille pas dans $\mathbb R$ mais sur un ensemble fini où on peut parler de successeur et de prédécesseur.
Nous travaillerons dans l'ensemble le plus commun, ${\mathbb V}_{64}$, celui des *binary64* c'est-à-dire des nombres en base 2 qui s'écrivent sur 64 bits selon le découpage suivant:

un bit de signe $s$ qui vaut 0 si le nombre est positif, 1 sinon ;
- 11 bits d'exposant décalé $e$ ;
- 53 bits de mantisse (ou « significande » en franglais) $ m, \quad 1 \leqslant m <2$.

Oups! $1 + 11 + 53 = 65$...En fait, on représente un nombre de ${\mathbb V}_{64}$ sous *forme normale* :

$$v = (-1)^s\times 1,f \times 2^E$$

avec $m=1,f$: point n'est besoin de stocker le 1 qui est implicite donc on ne stocke que $f$ et on gagne 1 bit. On appelle $f$ la pseudo-mantisse.

De même il faudrait, pour l'exposant, garder un bit pour stocker son signe. On va en fait le décaler pour que l'exposant stocké commence à 1 : on gagne encore un bit de stockage. Pour cela on le décale de l'exposant positif maximum - 1. On va voir que l'on réserve l'exposant 0 et l'exposant maximum pour représenter des nombres spéciaux.


En binary 64, cela signifie que le décalage vaut $\delta=11111111110_2=2^{11-1}-1=1023$ et $$e_{\rm stocké}=E_{\rm réel}+\delta$$


Pour en revenir à notre problème précédent, comment représenter $2^53+1$ sur 64 bits?
Sur machine $$2^{53} = (-1)^0\times 1,000...000 \times 2^{53+1023}$$

0 10000110100 0000000000000000000000000000000000000000000000000000

$$
1=(-1)^0\times 1,000...000 \times 2^{0+1023}
$$

0 01111111111 0000000000000000000000000000000000000000000000000000

Or, pour additionner deux nombres à virgule flottante, il faut d'abord exprimer le plus petit dans l'exposant du plus grand (donc aussi éventuellement le « dénormaliser »).

$$1=0,\underbrace{000...000}_{\rm 52 bits}\ 1 \times 2^{53}$$

Cela donne:




1,0000000000000000000000000000000000000000000000000000|


0,0000000000000000000000000000000000000000000000000000|1




1,0000000000000000000000000000000000000000000000000000|1




On ne peut pas mettre le dernier 1 dans la pseudo-mantisse qui n'a que 52 cases.

Le nombre obtenu n'est donc pas un « flottant » : il va falloir l'arrondir. La norme IEEE 754 propose quatre modes d'arrondi. Le mode par défaut est l'arrondi au plus proche (RN : Round to the Nearest). En cas d'équidistance, on choisit le flottant dont la pseudo-mantisse se termine par un 0.


Ici, il ne faut pas oublier que l'on travaille en base 2 : notre nombre, qui n'est pas un flottant de $\mathbb V_{64}$,
est à équidistance de $$1,\underbrace{000...000}{\rm 52 bits et } 1,\underbrace{000...001}{\rm 52 bits}$$
Il sera arrondi à celui qui se termine par un 0 : $2^{53}$.

$$
RN(2^{53}+1)=2^{53}
$$

Vous pouvez alors comprendre pourquoi:

$$
RN(2^{53}+3)=2^{53}+4
$$

La mauvaise nouvelle : une addition entre deux nombres flottants peut donc, dans certains cas, créer une erreur.


La bonne nouvelle : on peut récupérer cette erreur.

Les sommes compensées de deux flottants

Par la suite, on notera $x\ominus y$ l'arrondi de la somme de $x-y$ c'est-à-dire $x-y+err(x \ominus y)$.



Voici un lemme publié en 1973 par Pat Sterbenz :




Lemme de Sterbenz (1973)
Soit $(x,y)\in {\mathbb V}^2$ vérifiant $\dfrac{x}{2} \leqslant y \leqslant 2x$. On a
$$x\ominus y=x-y$$
La différence de deux nombres de $\mathbb V$ suffisamment proches est donc exacte.





Pour en avoir une démonstration, on se référera au poly pp 24 et 25.


La condition $\dfrac{x}{2} \leqslant y \leqslant 2x$ signifie en fait que $x$ et $y$ sont distants d'une puissance de 2 au maximum.
Que se passe-t-il quand on ne peut pas s'assurer à chaque instant que les
opérandes sont suffisamment proches?


Voici maintenant un premier algorithme dû à T.J. Dekker et W. Kahan qui répond à cette question:

Pseudo-langage
x,y : deux flottants tels que |x| >= |y|
s  <- x + y
yv <- s - x
d  <- y - yv
retourne (s,d) 

La ligne 1 donne $s=x + y + \operatorname{err}(x \oplus y)$.

De plus, comme $|x| \geqslant |y|$, c'est $x$ qui « gagne » donc $s$ et $x$
ont le même signe.

La ligne 2 calcule un $y$ virtuel ($y_v$) qui vaut $s \ominus x$: on aimerait
bien que ça fasse $s-x$ c'est-à-dire $y + \operatorname{err}(x \oplus y)$.

Enfin la ligne 3 calcule $y\ominus y_v$: on aimerait bien aussi que ça fasse
$y-y_v$ c'est-à-dire $- \operatorname{err}(x \oplus y)$ comme ça on récupérerait
exactement l'erreur commise.

En plus, comme $d=- \operatorname{err}(x \oplus y)=$, on a $|d| \leqslant \frac{1}{2}\operatorname{ulp}(s)$
donc $s$ et $d$ ne se chevauchent pas dans leur somme.

On aimerait bien, on aimerait bien...En fait tout s'arrange car nous allons
pouvoir utiliser le lemme de Sterbenz.

Il suffit de distinguer deux cas:

  • si $x$ et $y$ sont de même
    signe OU si $\mathbf{|y| \leqslant \frac{|x|}{2}}$: alors
    $\frac{x}{2}\leqslant s \leqslant 2x$ et on peut appliquer le lemme;
  • sinon: on a $x$ et $y$ de signes opposés ET $y>\frac{\vert x\vert}{2}$
    alors si $y$ est négatif $x/2 < -y < x$ et sinon -x/2<y<-x. Dans les deux
    cas, d'après le lemme de Sterbenz, $s$ est calculée exactement et
    alors $y_v=y$.

Les sommes compensées d'un nombre quelconque de flottants

L'idée de cet algorithme, dit de sommation en cascade, est d'additionner dans un accumulateur les erreurs puis les additionnée à la pseudo-somme à la fin. Cette idée a été publiée par Michèle Pichat en 1972 dans un mémoire de maîtrise et simultanément par Arnold Neumaier et ensuite généralisé en 2008 par Siegfried M. Rump, Takeshi Ogita
et Shin'Ichi Oishi.

L'idée tient en un petit dessin :

dessin1

On peut voir ce que cela donne sur un exemple:

dessin2

Par exemple avec Python

Python
def Fast2Sum(a,b):
    if a >= b:
        s = a + b
        z = s - a
        return (s, b - z)
    else:
        return Fast2Sum(b,a)
 
def Pichat(liste):
    accS = 0
    accE = 0
    for k in liste:
        s,e = Fast2Sum(accS,k)
        accS  = s
        accE += e
    return accS + accE

et on obtient:

Python
In [1]: Pichat([1.0 / k for k in range(1,100001)])
Out[1]: 12.090146129863427

Ou avec Haskell

Haskell
fastTwoSum a b  
    | abs a >= abs b = 
      let s = a + b in
      let z = s - a in 
      [s, b - z]
    | otherwise =  fastTwoSum b a
 
sommePichat liste =
  let [ss,ee] = foldl (\ [accS,accE] x -> 
                        let [s,e] = fastTwoSum accS x in [s,e + accE]) [0,0] liste 
  in  ss + ee

On peut comparer les trois sommes :

Haskell
*Main> sommePichat (map (\ k -> 1 / k) [1 .. 100000])
12.090146129863427
*Main> sum (map (\ k -> 1 / k) [1 .. 100000])
12.090146129863335
*Main> sum (map (\ k -> 1 / k) [100000, 99999 .. 1])
12.090146129863408

Sachant que le résultat correct est 12,0901461298634280...