/* $Id: pb-e.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 void calc_e (mpfr_t e, int prec)
{
  mpfr_t t, u;
  int i = 1;
#ifdef PREC
  int newprec;
#endif

  mpfr_inits2 (prec, t, u, (void *) 0);
  mpfr_set_ui (e, 2, MPFR_RNDN);
  mpfr_set_ui (t, 1, MPFR_RNDN);

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

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

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

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

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

  mpfr_init2 (e1, prec);
  mpfr_init2 (e2, prec + 8);
  mpfr_init2 (u, 2);

  c = clock ();
  calc_e (e1, prec);
  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: 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');

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

/* On prunille (2.7 GHz PowerPC G5, 32-bit ABI, Mac OS X 10.4.10):
   - without -DPREC: arg = 252000, error <= 3.47843e-75858
   - with -DPREC:    arg = 317600, error <= 1.65588e-95606
   On vin (3.0 GHz Intel Pentium D, Linux 2.6.22):
   - without -DPREC: arg = 228000, error <= 1.25554e-68633
   - with -DPREC:    arg = 261000, error <= 2.79019e-78567
   On mermoz (3.0 GHz Intel Pentium 4, Linux 2.6.22):
   - without -DPREC: arg = 171000, error <= 1.81990e-51475
   - with -DPREC:    arg = 194000, error <= 3.06676e-58398
   On courge (2.2 GHz AMD Opteron, Linux 2.6.22):
   - without -DPREC: arg = 331000, error <= 2.67530e-99640
   - with -DPREC:    arg = 374000, error <= 4.74827e-112584
   On zaurus (400 MHz Intel XScale-PXA255, Linux 2.4.18):
   - without -DPREC: arg = 67600, error <= 1.26836e-20348
   - with -DPREC:    arg = 78300, error <= 2.06511e-23569
*/