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

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

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

static void out (mpfr_ptr x, mpfr_ptr y)
{
  int i;

  for (i = 2; i <= 9; i++)
    {
      mpfr_pow_ui (y, x, i, MPFR_RNDN);
      printf ("1.001^%d = ", i);
      mpfr_out_str (stdout, 10, 28, y, MPFR_RNDN);
      putchar ('\n');
    }
}

int main (void)
{
  mpfr_t x, y;
  int i, j;

  mpfr_inits2 (100, x, y, (void *) 0);

  printf ("Incorrect results with mpfr_set_d:\n");
  mpfr_set_d (x, 1.001, MPFR_RNDN);
  out (x, y);

  printf ("\nWithout using double-precision arithmetic:\n");
  mpfr_set_ui (x, 1001, MPFR_RNDN);
  mpfr_div_ui (x, x, 1000, MPFR_RNDN);
  out (x, y);

  printf ("\nBounds:\n");
  for (i = 2; i <= 9; i++)
    {
      for (j = 0; j <= 1; j++)
        {
          printf (j ? " , " : "[ ");
          mpfr_set_ui (x, 1001, rnd[j]);
          mpfr_div_ui (x, x, 1000, rnd[j]);
          mpfr_pow_ui (y, x, i, rnd[j]);
          mpfr_out_str (stdout, 10, 28, y, rnd[j]);
        }
      printf (" ]\n");
    }

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

/*
Incorrect results with mpfr_set_d (on Linux/x86):
1.001^2 = 1.002000999999999779511483666
1.001^3 = 1.003003000999999668936492725
1.001^4 = 1.004006004000999558140572290
1.001^5 = 1.005010010005000447123391078
1.001^6 = 1.006015020015005336884617363
1.001^7 = 1.007021035035020231424918977
1.001^8 = 1.008028056070055140748964310
1.001^9 = 1.009036084126125084871426308

Without using double-precision arithmetic:
1.001^2 = 1.002001000000000000000000000
1.001^3 = 1.003003001000000000000000000
1.001^4 = 1.004006004001000000000000000
1.001^5 = 1.005010010005001000000000000
1.001^6 = 1.006015020015006001000000000
1.001^7 = 1.007021035035021007001000000
1.001^8 = 1.008028056070056028008001000
1.001^9 = 1.009036084126126084036009001

Bounds:
[ 1.002000999999999999999999999 , 1.002001000000000000000000001 ]
[ 1.003003000999999999999999999 , 1.003003001000000000000000001 ]
[ 1.004006004000999999999999999 , 1.004006004001000000000000001 ]
[ 1.005010010005000999999999999 , 1.005010010005001000000000001 ]
[ 1.006015020015006000999999999 , 1.006015020015006001000000001 ]
[ 1.007021035035021007000999999 , 1.007021035035021007001000001 ]
[ 1.008028056070056028008000999 , 1.008028056070056028008001001 ]
[ 1.009036084126126084036009000 , 1.009036084126126084036009002 ]
*/