/* $Id: pb-ex-sqrt.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, s, r;
  int iprec, oprec;

  if (argc != 5)
    {
      fprintf (stderr, "Usage: pb-ex-sqrt <x> <y> <iprec> <oprec>\n");
      exit (EXIT_FAILURE);
    }

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

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

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

  if (mpfr_set_str (x, argv[1], 0, MPFR_RNDN))
    {
      fprintf (stderr, "pb-ex-sqrt: bad value x (%s)\n", argv[1]);
      exit (EXIT_FAILURE);
    }

  if (mpfr_set_str (y, argv[2], 0, MPFR_RNDN))
    {
      fprintf (stderr, "pb-ex-sqrt: bad value y (%s)\n", argv[2]);
      exit (EXIT_FAILURE);
    }

  /* TODO: error analysis. Let's assume that computing x + y on
     oprec + 1 is sufficient. */
  mpfr_init2 (s, oprec + 1);
  mpfr_init2 (r, oprec);

  mpfr_clear_flags ();
  mpfr_add (s, x, y, MPFR_RNDN);

  if (mpfr_nan_p (s) || mpfr_cmp_ui (s, 0) < 0)
    mpfr_set_nan (r);
  else if (mpfr_overflow_p ())
    {
      mpfr_t t;

      printf ("x + y overflows...\n");
      if (mpfr_cmp (x, y) < 0)
        mpfr_swap (x, y);
      mpfr_sqrt (s, x, MPFR_RNDN);
      mpfr_init2 (t, oprec + 2);
      mpfr_div (t, y, x, MPFR_RNDN);
      mpfr_add_ui (t, t, 1, MPFR_RNDN);
      mpfr_sqrt (t, t, MPFR_RNDN);
      mpfr_mul (r, s, t, MPFR_RNDN);
      mpfr_clear (t);
    }
  else if (mpfr_underflow_p ())
    {
      int logscale;

      printf ("x + y underflows...\n");
      /* Unless the precisions are huge, we can assume that the following
         scaling by 2^iprec or 2^(iprec-1) won't lead to an overflow. */
      logscale = iprec >> 1;
      mpfr_mul_2ui (x, x, logscale << 1, MPFR_RNDN);
      mpfr_mul_2ui (y, y, logscale << 1, MPFR_RNDN);
      /* ulp(x) and ulp(y) are now representable. */
      mpfr_add (s, x, y, MPFR_RNDN);  /* no possible underflow */
      mpfr_sqrt (r, s, MPFR_RNDN);
      mpfr_div_2ui (r, r, logscale, MPFR_RNDN);
    }
  else
    mpfr_sqrt (r, s, MPFR_RNDN);

  mpfr_out_str (stdout, 16, 0, r, MPFR_RNDN);
  putchar ('\n');

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

/*
$ ./pb-ex-sqrt 144 145 32 32
1.10000000@1
$ ./pb-ex-sqrt 144 145 32 128
1.10000000000000000000000000000000@1

$ ./pb-ex-sqrt 4611686018427387904 4294967298 32 32
8.00000010@7
$ ./pb-ex-sqrt 4611686018427387904 4294967298 32 128
8.000000100000000fffffffe000000030@7

$ ./pb-ex-sqrt 0x.7fffffffe@268435456 17 34 32
x + y overflows...
b.504f3340@134217727
$ ./pb-ex-sqrt 0x.7fffffffe@268435456 17 34 128
b.504f333e33dc61dd8d7b33c81fa0f630@134217727

$ ./pb-ex-sqrt 0x.7fffffffe@268435456 0x.1@268435456 34 32
x + y overflows...
c.00000000@134217727
$ ./pb-ex-sqrt 0x.7fffffffe@268435456 0x.1@268435456 34 128
x + y overflows...
b.fffffffeaaaaaaaa97b425ed075fde50@134217727

$ ./pb-ex-sqrt 0x.1ffffffff8@-268435455 -0x.1ffff@-268435455 34 32
x + y underflows...
f.fffc0000@-134217731
$ ./pb-ex-sqrt 0x.1ffffffff8@-268435455 -0x.1ffff@-268435455 34 128
x + y underflows...
f.fffbffff7fffdffff5fffc7ffeafff80@-134217731

$ ./pb-ex-sqrt @Inf@ 1 2 32
@Inf@
*/