/*
 * $Id: rump.c 65304 2013-11-28 12:59:24Z vinc17/ypig $
 *
 * Link with -lmpfr -lgmp -lm *in that order*
 *
 * Computes
 *   333.75 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) + 5.5 b^8 + a / (2b)
 * with a = 77617 and b = 33096 ("Rump's polynomial").
 *
 * Same semantics as in C with correct rounding,
 * but with different precisions.
 *
 * Arguments: mode, precmin and precmax.
 * Display mode 0: full number of digits.
 * Display mode 1: more compact form.
 */

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

static mpfr_prec_t readprec (const char *str, const char *var)
{
  long int prec;
  char *end;

  prec = strtol (str, &end, 10);

  if (*end != '\0')
    {
      fprintf (stderr, "rump: error when reading %s\n", var);
      exit (EXIT_FAILURE);
    }

  if (prec < MPFR_PREC_MIN)
    {
      fprintf (stderr, "rump: %s is too small\n", var);
      exit (EXIT_FAILURE);
    }

  if (prec > MPFR_PREC_MAX || errno == ERANGE)
    {
      fprintf (stderr, "rump: %s is too large\n", var);
      exit (EXIT_FAILURE);
    }

  return prec;
}

static const mpfr_rnd_t rnd[3] = { MPFR_RNDD, MPFR_RNDN, MPFR_RNDU };
static const char *rndname[3] =
  { "Lower bound", "Using MPFR_RNDN", "Upper bound" };
static int mode;

static void rump (mpfr_prec_t prec)
{
  mpfr_t x[24], t1, t2;
  mpfr_t *const a  = x + 3;
  mpfr_t *const b  = x + 6;
  mpfr_t *const a2 = x + 9;
  mpfr_t *const b2 = x + 12;
  mpfr_t *const b4 = x + 15;
  mpfr_t *const b6 = x + 18;
  mpfr_t *const b8 = x + 21;
  int i, j;

  mpfr_set_default_prec (prec);
  switch (mode)
    {
    case 0:
      printf ("Precision %ld:\n", prec);
      break;
    case 1:
      printf ("%8ld   ", prec);
      break;
    }

  for (i = 0; i < sizeof(x) / sizeof(mpfr_t); i++)
    mpfr_init (x[i]);
  mpfr_init (t1);
  mpfr_init (t2);

  for (i = 0; i < 3; i++)
    {
      mpfr_set_ui (a[i], 77617, rnd[i]);
      mpfr_set_ui (b[i], 33096, rnd[i]);
      mpfr_pow_ui (a2[i], a[i], 2, rnd[i]);
      for (j = 1; j <= 4; j++)
        mpfr_pow_ui (b2[3*(j-1)+i], b[i], 2*j, rnd[i]);
    }

  for (i = 0; i < 3; i++)
    {
      mpfr_set_d (t1, 333.75, rnd[i]);          // t1: 333.75
      mpfr_mul (x[i], t1, b6[i], rnd[i]);       // x:  333.75 b^6
      mpfr_mul_ui (t1, a2[i], 11, rnd[i]);      // t1: 11 a^2
      mpfr_mul (t1, t1, b2[i], rnd[i]);         // t1: 11 a^2 b^2
      mpfr_sub (t1, t1, b6[2-i], rnd[i]);       // t1: 11 a^2 b^2 - b^6
      mpfr_mul_ui (t2, b4[2-i], 121, rnd[2-i]); // t2: 121 b^4
      mpfr_sub (t1, t1, t2, rnd[i]);            // t1: 11 a^2 ... - 121 b^4
      mpfr_sub_ui (t1, t1, 2, rnd[i]);          // t1: 11 a^2 b^2 - ... - 2
      mpfr_mul (t1, t1, a2[MPFR_SIGN(t1) < 0 ? 2-i : i], rnd[i]);
      mpfr_add (x[i], x[i], t1, rnd[i]);        // x:  333.75 b^6 + a^2 (...)
      mpfr_set_d (t1, 5.5, rnd[i]);             // t1: 5.5
      mpfr_mul (t1, t1, b8[i], rnd[i]);         // t1: 5.5 b^8
      mpfr_add (x[i], x[i], t1, rnd[i]);        // x: 333.75 ... + 5.5 b^8
      mpfr_mul_2ui (t1, b[2-i], 1, rnd[2-i]);   // t1: 2b
      mpfr_div (t1, a[i], t1, rnd[i]);          // t1: a / (2b)
      mpfr_add (x[i], x[i], t1, rnd[i]);        // x: 333.75 ... + a / (2b)

      switch (mode)
        {
        case 0:
          printf ("%16s: ", rndname[i]);
          mpfr_out_str (stdout, 10, 0, x[i], rnd[i]);
          putchar ('\n');
          break;
        case 1:
          {
            size_t len;
            len = mpfr_out_str (stdout, 10, 15, x[i], rnd[i]);
            if (i < 2)
              while (len++ < 23)
                putchar (' ');
          }
          break;
        }
    }

  if (mode == 1)
    putchar ('\n');

  for (i = 0; i < sizeof(x) / sizeof(mpfr_t); i++)
    mpfr_clear (x[i]);
  mpfr_clear (t1);
  mpfr_clear (t2);
}

int main (int argc, char *argv[])
{
  mpfr_prec_t precmin, precmax;

  if (argc != 4)
    {
      fprintf (stderr, "Usage: rump <mode> <precmin> <precmax>\n");
      exit (EXIT_FAILURE);
    }

  mode = atoi (argv[1]);
  if (mode != 0 && mode != 1)
    {
      fprintf (stderr, "rump: mode must be 0 or 1\n");
      exit (EXIT_FAILURE);
    }

  precmin = readprec (argv[2], "precmin");
  precmax = readprec (argv[3], "precmax");

  if (mode == 1)
    printf ("Precision, lower bound, value using MPFR_RNDN, upper bound:\n");

  while (precmin <= precmax)
    rump (precmin++);

  return 0;
}