Problème — Sommer un million de termes, et l'ordre qui change tout

1. Les deux sommes, mesurées.

n = 10 000 000
float, k croissant   : 15,403682709
float, k decroissant : 16,686031342
double (reference)   : 16,695311366
ln(n) + gamma        : 16,695311316      (valeur theorique)

erreur, ordre croissant   : 1,291628657
erreur, ordre decroissant : 0,009280024

Le même calcul, les mêmes $10^7$ additions, les mêmes termes : deux résultats qui diffèrent de $1,28$, soit près de $8\%$. Et l'ordre décroissant est $139$ fois plus précis que l'ordre croissant.

La colonne de référence n'est pas une commodité : la somme en double coïncide avec $\ln n+\gamma$ sur sept décimales, où $\gamma \approx 0,5772156649$ est la constante d'Euler. C'est ce qui autorise à parler d'« erreur » pour les deux autres.

2. L'explication : l'absorption. Dans l'ordre croissant, la somme partielle atteint vite une taille de l'ordre de $10$, tandis que les termes $1/k$ deviennent minuscules. Or un float ne porte que $24$ bits de mantisse : à côté d'un total $s\approx 15\approx 2^{3,9}$, l'écart entre deux flottants consécutifs vaut environ $2^4\times 2^{-24}=2^{-20}\approx 10^{-6}$. Tout terme inférieur à la moitié de cet écart est purement et simplement absorbé : $s+1/k$ rend $s$.

Dans l'ordre décroissant, les petits termes sont additionnés entre eux d'abord, quand le total est encore petit. Ils s'accumulent au lieu de disparaître, et ne rejoignent le gros total qu'une fois devenus significatifs.

3. Le rang d'absorption, mesuré :

premier k tel que s + 1/k == s (en float) : 2097152

Et $2097152=2^{21}$, exactement la valeur prédite par le raisonnement ci-dessus ($1/k<2^{-21}$ équivaut à $k>2^{21}$). Sur les $10^7$ termes, près de $8$ millions n'ont donc contribué à rien du tout dans la version croissante — le programme a passé la majeure partie de son temps à ajouter zéro.

4. La règle. Quand on somme des termes de tailles très différentes, on commence par les plus petits. Plus généralement :

• on n'additionne pas un petit nombre à un grand si on peut l'éviter ;
• on préfère un double à un float dès qu'on somme beaucoup de termes — la même expérience en double donne les deux ordres identiques sur neuf décimales ;
• et l'on garde à l'esprit que deux programmes mathématiquement équivalents ne le sont pas numériquement. C'est la conséquence de la perte d'associativité mesurée à l'exercice 3.6.