HPC Magazine novembre 2013 - Des outils pour calculer avec précision.

Listing 1 - Calcul de la formule de Rump en C classique.

#include <math.h>
#include <stdio.h>

main()
{
  float a = 77617.;
  float b = 33096.;
  float res;

  res = (333.0 + 3.0/4.0)*(b*b*b*b*b*b) + (a*a)*(11.0*(a*a)*(b*b) - (b*b*b*b*b*b) - 
    121.0*(b*b*b*b) -2.0) + 11.0/2.0*(b*b*b*b*b*b*b*b) + a/(2.0*b);
  printf("res=%.14e\n",res);

  res=333.75*b*b*b*b*b*b+a*a*(11.0*a*a*b*b-b*b*b*b*b*b-121.0*b*b*b*b-2.0) +
    5.5*b*b*b*b*b*b*b*b+a/(2.0*b);
  printf("res=%.14e\n",res);

  /* On my x86 64-bit AMD, the results are:
  /* in single precision: res=-2.34119289742886e+29, res=7.08931176699221e+29
  /* in double precision : res=-1.18059162071741e+21, res=1.17260394005318e+00  */
}