/* Ecole CNC'2
   Nancy 25-26 juin 2009

   probleme 1.1 */

/* include stdio.h before including mpfr.h when mpfr_printf or mpfr_out_str is
   needed */
#include <stdio.h>
#include <mpfr.h>

void
sum (mpfr_srcptr x)
{
  unsigned int i;
  mpfr_t s, t, u;
  mp_prec_t prec = mpfr_get_prec (x);

  mpfr_init2 (t, prec);
  mpfr_set_ui (t, 1, GMP_RNDD);
  mpfr_init2 (s, prec);
  mpfr_set_ui (s, 1, GMP_RNDD);
  mpfr_init2 (u, prec);
  for (i = 1; i <= 100; i++)
    {
      mpfr_mul_ui (t, t, i, GMP_RNDU);
      mpfr_pow_ui (u, x, i, GMP_RNDD);
      mpfr_div (u, u, t, GMP_RNDD);
      mpfr_add (s, s, u, GMP_RNDD);
    }
  printf ("Sum is ");
  mpfr_out_str (stdout, 10, 0, s, GMP_RNDN);
  putchar ('\n');

  mpfr_clear (s);
  mpfr_clear (t);
  mpfr_clear (u);
  return;
}

int
main (void)
{
  mpfr_t x, y;
  mp_prec_t p;

  mpfr_init (y);

  /* x := 1.5 */
  mpfr_init (x);
  mpfr_set_d (x, 1.5, GMP_RNDN);
  
  mpfr_set_prec (y, 150);

  mpfr_printf ("y = 1 + 1.5^1/1! + 1.5^2/2! + ... + 1.5^100/100!\n       ");
  sum (x);

  mpfr_printf ("y = mpfr_exp(1.5)\n");

  mpfr_exp (y, x, GMP_RNDN);
  mpfr_printf ("y[150 bits] = %RNf\n", y);

  mpfr_set_prec (y, 170);
  mpfr_exp (y, x, GMP_RNDN);
  mpfr_printf ("y[170 bits] = %RNf\n", y);

  printf ("---\n");

  for (p = 130; p < 151; p++)
    {
      mpfr_set_prec (y, p);
      mpfr_exp (y, x, GMP_RNDN);

      mpfr_printf ("y[%Pu bits] = %.41RNf\n", p, y);
    }

  mpfr_clears (x, y, (mpfr_t*)0);
 
  return 0;
}