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

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

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

  if (argc != 5 && argc != 6)
    {
      fprintf (stderr, "Usage: pb-cr-canround <x> <xprec> <n> <yprec>"
               " [ <no_inex> ]\n");
      exit (EXIT_FAILURE);
    }
  reqinex = argc == 5;

  xprec = atoi (argv[2]);
  if (xprec < MPFR_PREC_MIN || xprec > MPFR_PREC_MAX)
    {
      fprintf (stderr, "pb-cr-canround: 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-canround: bad value x (%s)\n", argv[1]);
      exit (EXIT_FAILURE);
    }

  n = atoi (argv[3]);
  if (n < 2)
    {
      fprintf (stderr, "pb-cr-canround: 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-canround: incorrect output precision %d\n",
               yprec);
      exit (EXIT_FAILURE);
    }

  for (rnd = 0; rnd < 4; rnd++)
    {
      int inex;
      mp_prec_t tprec = yprec;

      mpfr_init2 (y, tprec);
      do
        {
          tprec += 8;
          mpfr_set_prec (y, tprec);
          printf ("Intermediate precision %ld...\n", (long) tprec);
          mpfr_clear_inexflag ();
          mpfr_root (y, x, n, MPFR_RNDN);
          mpfr_ui_div (y, 1, y, MPFR_RNDN);
        }
      while (mpfr_inexflag_p () && ! mpfr_can_round (y, tprec - 2, MPFR_RNDN,
          reqinex ? MPFR_RNDZ : rnd, yprec + (reqinex && rnd == MPFR_RNDN)));
      inex = mpfr_prec_round (y, yprec, rnd);
      printf ("%s: ", mpfr_print_rnd_mode (rnd));
      mpfr_out_str (stdout, 2, 0, y, rnd);
      printf (reqinex ? " (inex = %d)\n" : "\n", inex < 0 ? -1 : inex > 0);
      mpfr_clear (y);
    }

  mpfr_clear (x);
  return 0;
}

/*
Let t = root(x,n), y1 = o(t), u = 1/t, y2 = o(1/y1).
y1 = t(1 + Eps), where Eps = [-eps,eps] and eps = 2^(-p).
1/y1 = y2(1 + Eps), thus u = 1/t = y2(1 + Eps)^2.
|y2 - u| <= y2 (2 eps + eps^2) <= 2^(E(y2) + 2 - p).

$ ./pb-cr-canround 12 32 17 31 -
Intermediate precision 39...
Intermediate precision 47...
MPFR_RNDN: 1.101110100101111110000011100011e-1
Intermediate precision 39...
MPFR_RNDZ: 1.101110100101111110000011100010e-1
Intermediate precision 39...
MPFR_RNDU: 1.101110100101111110000011100011e-1
Intermediate precision 39...
MPFR_RNDD: 1.101110100101111110000011100010e-1
$ ./pb-cr-canround 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-canround 12 32 17 32 -
Intermediate precision 40...
MPFR_RNDN: 1.1011101001011111100000111000101e-1
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDZ: 1.1011101001011111100000111000101e-1
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDU: 1.1011101001011111100000111000110e-1
Intermediate precision 40...
Intermediate precision 48...
MPFR_RNDD: 1.1011101001011111100000111000101e-1
$ ./pb-cr-canround 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)
*/