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

Listing 3 - Calcul de la formule de Rump avec la bibliothèque MPFI.

#include <mpfi.h>
#include <mpfi_io.h>

int main(int argc, char *argv[])
{
  int prec=53;
  int go_on;
  char ch;
  mpfi_t a, b, c, d, e, f, g, h, i;

  mpfr_set_default_prec(prec);
  mpfi_init(a); mpfi_init(b); mpfi_init(c);
  mpfi_init(d); mpfi_init(e); mpfi_init(f);
  mpfi_init(g); mpfi_init(h); mpfi_init(i);

  do {
    /* Setting the desired precision */
    printf("Enter the computing precision: ");
    scanf("%d", &prec);
    printf("\t%d\n", prec);
    mpfr_set_default_prec(prec);
    mpfi_set_prec(a,prec); mpfi_set_prec(b,prec); mpfi_set_prec(c,prec);
    mpfi_set_prec(d,prec); mpfi_set_prec(e,prec); mpfi_set_prec(f,prec);
    mpfi_set_prec(g,prec); mpfi_set_prec(h,prec); mpfi_set_prec(i,prec);

    /* Computing the different terms of Rump’s formula */
    mpfi_set_si(a, 77617);
    mpfi_set_si(b, 33096);
    mpfi_set_si(c, 3); mpfi_div_si(c,c,4); mpfi_add_si(c,c,333);
    mpfi_sqr(i,b);            /* i = b^2 */
    mpfi_sqr(h,i);            /* h = b^4 */
    mpfi_mul(g, i, h);        /* g = b^6 */
    mpfi_mul(c, c, g);        /* Here c = (333+3/4)*b^6 */

    mpfi_sqr(f, a);           /* f = a^2 */
    mpfi_mul_si(d, f, 11);    /* d = 11*a^2 */
    mpfi_mul(d, d, i);        /* d = 11*a^2*b^2 */
    mpfi_sub(d, d, g);        /* d = 11*a^2*b^2 -b^6 */
    mpfi_mul_si(e, h, 121);   /* e = 121*b^4 */
    mpfi_sub(d, d, e);        /* d = 11*a^2*b^2 -b^6 -121*b^4 */
    mpfi_sub_si(d, d, 2);     /* d = 11*a^2*b^2 -b^6 -121*b^4 -2 */
    mpfi_mul(d, d, f);        /* d = (11*a^2*b^2 -b^6 -121*b^4 -2)*a^2 */
    mpfi_add(c, c, d);
    mpfi_sqr(f, a);           /* f = a^2 */
    mpfi_mul_si(d, f, 11);    /* d = 11*a^2 */
    mpfi_mul(d, d, i);        /* d = 11*a^2*b^2 */
    mpfi_sub(d, d, g);        /* d = 11*a^2*b^2 -b^6 */
    mpfi_mul_si(e, h, 121);   /* e = 121*b^4 */
    mpfi_sub(d, d, e);        /* d = 11*a^2*b^2 -b^6 -121*b^4 */
    mpfi_sub_si(d, d, 2);     /* d = 11*a^2*b^2 -b^6 -121*b^4 -2 */
    mpfi_mul(d, d, f);        /* d = (11*a^2*b^2 -b^6 -121*b^4 -2)*a^2 */
    mpfi_add(c, c, d);

    mpfi_sqr(h, h);           /* h = b^8 */
    mpfi_mul_si(e, h, 11);    /* e = 11*b^8 */
    mpfi_div_si(e, e, 2);     /* e = 11/2*b^8 */
    mpfi_add(c, c, e);

    mpfi_div(f, a, b);        /* f = a/b */
    mpfi_div_si(f, f, 2);     /* f = a/(2b) */
    mpfi_add(c, c, f);

    /* Printing results */
    printf("\nEvaluation =  ");
    mpfi_out_str(stdout, 10, 0, c);

    /* Try again? */
    printf("\nDo you want more Rump’s formula? y/n ");
    do {scanf("%c", &ch);} while (isspace(ch));
    go_on=(ch==’y’); printf("\n");
  } while (go_on) ;

  mpfi_clear(a); mpfi_clear(b); mpfi_clear(c);
  mpfi_clear(d); mpfi_clear(e); mpfi_clear(f);
  mpfi_clear(g); mpfi_clear(h); mpfi_clear(i);
}