/*
 * $Id: pb-muller2.c 47490 2011-11-10 17:07:08Z vinc17/ypig $
 *
 * Link with -lmpfr -lgmp *in that order*
 *
 * Example by Jean-Michel Muller:
 *   u_0 = e - 1
 *   u_n = n u_{n-1} - 1
 *
 * Arguments: precision, n
 */

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>

static const mp_rnd_t rnd[3] = { MPFR_RNDD, MPFR_RNDN, MPFR_RNDU };

static void display (mpfr_t x[3], long k)
{
  int i;

  printf ("%6ld   ", k);
  for (i = 0; i < 3; i++)
    {
      size_t len;
      len = mpfr_out_str (stdout, 10, 17, x[i], rnd[i]);
      if (i < 2)
        while (len++ < 24)
          putchar (' ');
    }
  putchar ('\n');
}

int main (int argc, char *argv[])
{
  int prec, n, i, k, ndec = 0;
  mpfr_t u[3], one;

  if (argc != 3 && argc != 4)
    {
      fprintf (stderr, "Usage: pb-muller2 <prec> <n> [ <ndec> ]\n");
      exit (EXIT_FAILURE);
    }

  prec = atoi (argv[1]);
  if (prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX)
    {
      fprintf (stderr, "pb-muller2: incorrect precision %d\n", prec);
      exit (EXIT_FAILURE);
    }
  mpfr_set_default_prec (prec);

  n = atoi (argv[2]);
  if (n <= 0)
    {
      fprintf (stderr, "pb-muller2: n must be positive\n");
      exit (EXIT_FAILURE);
    }

  if (argc > 3)
    {
      ndec = atoi (argv[3]);
      if (ndec < 2)
        {
          fprintf (stderr, "pb-muller2: ndec must be at least 2\n");
          exit (EXIT_FAILURE);
        }
    }

  mpfr_init2 (one, MPFR_PREC_MIN);
  mpfr_set_ui (one, 1, MPFR_RNDN);

  for (i = 0; i < 3; i++)
    {
      mpfr_init (u[i]);
      mpfr_expm1 (u[i], one, rnd[i]);
    }

  if (!ndec)
    display (u, 0);

  for (k = 1; k <= n; k++)
    {
      for (i = 0; i < 3; i++)
        {
          mpfr_mul_ui (u[i], u[i], k, rnd[i]);
          mpfr_sub_ui (u[i], u[i], 1, rnd[i]);
        }
      if (!ndec)
        display (u, k);
    }

  if (ndec)
    {
      printf ("Lower bound, value using MPFR_RNDN, upper bound:\n");
      for (i = 0; i < 3; i++)
        {
          mpfr_out_str (stdout, 10, ndec, u[i], rnd[i]);
          putchar ('\n');
        }
    }

  for (i = 0; i < 3; i++)
    mpfr_clear (u[i]);
  mpfr_clear (one);

  return 0;
}

/*
$ ./pb-muller2 8586 998 30
Lower bound, value using MPFR_RNDN, upper bound:
1.00200300300200300498581754171e-3
1.00200300300200300500412238262e-3
1.00200300300200300500412238263e-3
^^^^^^^^^^^^^^^^^^^^^^^ first 22 decimals
$ ./pb-muller2 8587 998 30
Lower bound, value using MPFR_RNDN, upper bound:
1.00200300300200300499496996216e-3
1.00200300300200300499496996217e-3
1.00200300300200300500412238263e-3
^^^^^^^^^^^^^^^^^^^^ first 19 decimals
$ ./pb-muller2 8588 998 30
Lower bound, value using MPFR_RNDN, upper bound:
1.00200300300200300499496996216e-3
1.00200300300200300499496996217e-3
1.00200300300200300499954617240e-3
^^^^^^^^^^^^^^^^^^^^ first 19 decimals
The answer is 1.002003003002003004e-3 with prec >= 8587, and >= 8588 to
guarantee the result in interval arithmetic, if mpfr_expm1 is used.
*/