/* $Id: pb-e-dir.c 47490 2011-11-10 17:07:08Z vinc17/ypig $ */

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

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

static void calc_e (mpfr_t e, int prec, int r)
{
  mpfr_t t, u;
  int i = 1;
#ifdef PREC
  int newprec;
#endif

  mpfr_inits2 (prec, t, u, (void *) 0);
  mpfr_set_ui (e, 2, rnd[r]);
  mpfr_set_ui (t, 1, rnd[r]);

  while (2 - mpfr_get_exp (t) < prec)
    {
      mpfr_div_ui (t, t, ++i, rnd[r]);
#ifdef PREC
      /* make -W pb-e-dir.c pb-e-dir CPPFLAGS=-DPREC */
      newprec = prec + mpfr_get_exp (t);
      if (newprec < MPFR_PREC_MIN)
        break;
      mpfr_prec_round (t, prec + mpfr_get_exp (t), rnd[r]);
#endif
      mpfr_add (e, e, t, rnd[r]);
    }

  mpfr_clears (t, u, (void *) 0);
}

int main (int argc, char *argv[])
{
  mpfr_t e1, e2, u, einf, esup, length;
  int prec, inex;
  clock_t c;
  double t;

  if (argc != 2)
    {
      fprintf (stderr, "Usage: pb-e-dir <prec>\n");
      exit (EXIT_FAILURE);
    }

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

  mpfr_inits2 (prec, e1, einf, esup, (void *) 0);
  mpfr_init2 (e2, prec + 8);
  mpfr_init2 (u, 2);

  c = clock ();
  calc_e (e1, prec, 1);
  t = (double) (clock () - c) / CLOCKS_PER_SEC;
  printf ("Time: %.2f s\n", t);

  mpfr_set_ui (u, 1, MPFR_RNDN);
  mpfr_exp (e2, u, MPFR_RNDN);

  mpfr_set_ui_2exp (u, 1, mpfr_get_exp (e2) - (mp_exp_t) mpfr_get_prec (e2)
                    - 1, MPFR_RNDN);  /* 1/2 ulp(e2), error bound for e2 */
  inex = mpfr_sub (e2, e2, e1, MPFR_RNDZ);  /* should be exact as e1 ~= e2 */
  mpfr_abs (e2, e2, MPFR_RNDN);
  if (inex)
    {
      fprintf (stderr, "pb-e-dir: warning: e1 -e2 is inexact\n");
      mpfr_nextabove (e2);  /* upper bound on initial |e1 - e2| */
    }
  mpfr_add (e2, e2, u, MPFR_RNDU);  /* upper bound on the error |e1 - e| */
  printf ("Error           <= ");
  mpfr_out_str (stdout, 10, 6, e2, MPFR_RNDU);
  putchar ('\n');

  calc_e (einf, prec, 0);
  calc_e (esup, prec, 2);
  mpfr_init2 (length, 32);
  inex = mpfr_sub (length, esup, einf, MPFR_RNDU);
  if (inex)
    fprintf (stderr, "pb-e-dir: warning: esup -einf is inexact\n");
  printf ("Interval length <= ");
  mpfr_out_str (stdout, 10, 6, length, MPFR_RNDU);
  putchar ('\n');

  mpfr_clears (e1, e2, u, einf, esup, length, (void *) 0);
  return 0;
}