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

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

#define SIGN(I) ((I) < 0 ? -1 : (I) > 0)

int main (int argc, char *argv[])
{
  mpfr_t x, y, z, t;
  int n, xprec, yprec, rnd;

  if (argc != 5)
    {
      fprintf (stderr, "Usage: pb-cr-interv <x> <xprec> <n> <yprec>\n");
      exit (EXIT_FAILURE);
    }

  xprec = atoi (argv[2]);
  if (xprec < MPFR_PREC_MIN || xprec > MPFR_PREC_MAX)
    {
      fprintf (stderr, "pb-cr-interv: incorrect input precision %d\n",
               xprec);
      exit (EXIT_FAILURE);
    }

  mpfr_init2 (x, xprec);
  if (mpfr_set_str (x, argv[1], 0, MPFR_RNDN) || mpfr_cmp_ui (x, 0) < 0)
    {
      fprintf (stderr, "pb-cr-interv: bad value x (%s)\n", argv[1]);
      exit (EXIT_FAILURE);
    }

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

  yprec = atoi (argv[4]);
  if (yprec < MPFR_PREC_MIN || yprec > MPFR_PREC_MAX)
    {
      fprintf (stderr, "pb-cr-interv: incorrect output precision %d\n",
               yprec);
      exit (EXIT_FAILURE);
    }

  for (rnd = 0; rnd < 4; rnd++)
    {
      int inex1, inex2;

      mpfr_inits2 (yprec, y, z, t, (void *) 0);
      do
        {
          mpfr_set_prec (t, mpfr_get_prec (t) + 8);
          printf ("Intermediate precision %d...\n", (int) mpfr_get_prec (t));
          inex2 = mpfr_root (t, x, n, MPFR_RNDZ);
          inex1 = mpfr_ui_div (y, 1, t, rnd);
          inex1 = SIGN (inex1);
          if (!inex1 && !inex2)
            break;  /* exact, i.e. power of two */
          mpfr_nextabove (t);
          inex2 = mpfr_ui_div (z, 1, t, rnd);
        }
      while (inex1 != SIGN (inex2) || ! mpfr_equal_p (y, z));
      printf ("%s: ", mpfr_print_rnd_mode (rnd));
      mpfr_out_str (stdout, 2, 0, y, rnd);
      printf (" (inex = %d)\n", inex1);
      mpfr_clears (y, z, t, (void *) 0);
    }

  mpfr_clear (x);
  return 0;
}

/*
$ ./pb-cr-interv 12 32 17 31
Intermediate precision 39...
Intermediate precision 47...
MPFR_RNDN: 1.101110100101111110000011100011e-1 (inex = 1)
Intermediate precision 39...
MPFR_RNDZ: 1.101110100101111110000011100010e-1 (inex = -1)
Intermediate precision 39...
MPFR_RNDU: 1.101110100101111110000011100011e-1 (inex = 1)
Intermediate precision 39...
MPFR_RNDD: 1.101110100101111110000011100010e-1 (inex = -1)
$ ./pb-cr-interv 12 32 17 32
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDN: 1.1011101001011111100000111000101e-1 (inex = -1)
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDZ: 1.1011101001011111100000111000101e-1 (inex = -1)
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDU: 1.1011101001011111100000111000110e-1 (inex = 1)
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDD: 1.1011101001011111100000111000101e-1 (inex = -1)
*/