/*
 * $Id: pb-rump.c 47491 2011-11-10 17:09:56Z vinc17/ypig $
 *
 * Link with -lmpfr -lgmp *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.
 * Same semantics as in C with exact rounding, but with different precisions.
 * Arguments: precmin and precmax.
 */

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

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

static void rump (mp_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);
  printf ("%8ld   ", prec);

  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++)
    {
      size_t len;

      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)

      len = mpfr_out_str (stdout, 10, 15, x[i], rnd[i]);
      if (i < 2)
        while (len++ < 23)
          putchar (' ');
    }

  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[])
{
  int precmin, precmax;

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

  precmin = atoi (argv[1]);
  precmax = atoi (argv[2]);
  if (precmin < MPFR_PREC_MIN || precmax < precmin ||
      precmax > MPFR_PREC_MAX)
    {
      fprintf (stderr, "pb-rump: incorrect precisions %d and %d\n",
               precmin, precmax);
      exit (EXIT_FAILURE);
    }

  printf ("Precision      Lower bound        "
          "Value using MPFR_RNDN      Upper bound\n");

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

  return 0;
}

/*
Precision      Lower bound        Value using MPFR_RNDN      Upper bound
      17   -4.86777830487641e32   1.17260742187500       4.05653143833191e32
      18   -2.43388915243821e32   1.17260742187500       2.02825333976556e32
      19   -1.21694457621911e32   1.17260360717773       1.01412357503269e32
      20   -6.08472288109551e31   1.17260360717773       4.05648965785558e31
      21   -3.04236144054776e31   1.17260360717773       2.53530313473778e31
      22   -1.52118072027388e31   1.17260408401489       1.26765108379856e31
      23   -7.60590360136938e30   1.17260384559631       7.60590481029520e30
      24   -3.16912650057058e30   -6.33825300114115e29   3.16912680280203e30
      25   -1.90147590034235e30   1.17260396480560       1.58456332584316e30
      26   -1.10919427519971e30   -1.58456325028529e29   6.33825319003581e29
      27   -4.75368975085587e29   1.17260393500328       3.16912654779424e29
      28   -2.37684487542794e29   -7.92281625142643e28   1.98070407466253e29
      29   -1.58456325028529e29   1.98070406285661e28    7.92281628094123e28
      30   -5.94211218856983e28   -9.90352031428304e27   4.95176016452022e28
      31   -3.46623210999907e28   4.95176015714152e27    1.98070406470129e28
      32   -1.48552804714246e28   2.47588007857076e27    1.48552804760363e28
      33   -6.18970019642691e27   -1.23794003928538e27   6.18970019757983e27
      34   -3.71382011785615e27   1.23794003928538e27    2.47588007885900e27
      35   -1.23794003928539e27   -3.09485009821345e26   1.85691005900013e27
      36   -7.73712524553363e26   1.17260394006735       7.73712524571378e26
      37   -3.86856262276682e26   -1.54742504910673e26   3.86856262281185e26
      38   -2.32113757366009e26   7.73712524553363e25    1.16056878683568e26
      39   -9.67140655691704e25   -1.93428131138341e25   9.67140655694519e25
      40   -4.83570327845852e25   9.67140655691703e24    3.86856262277386e25
      41   -2.41785163922926e25   1.17260394005280       1.93428131138517e25
      42   -1.20892581961463e25   1.17260394005325       1.20892581961507e25
      43   -7.25355491768778e24   1.17260394005325       6.04462909807425e24
      44   -3.02231454903658e24   -6.04462909807315e23   2.41785163922954e24
      45   -1.81338872942195e24   -3.02231454903657e23   1.20892581961470e24
      46   -1.05781009216281e24   1.17260394005316       4.53347182355495e23
      47   -6.04462909807315e23   7.55578637259143e22    2.26673591177746e23
      48   -2.26673591177743e23   1.17260394005318       1.51115727451830e23
      49   -1.13336795588872e23   1.88894659314786e22    7.55578637259146e22
      50   -4.72236648286965e22   -9.44473296573929e21   4.72236648286966e22
      51   -2.83341988972179e22   1.17260394005318       1.88894659314786e22
      52   -1.41670994486090e22   1.17260394005318       9.44473296573930e21
      53   -5.90295810358706e21   -1.18059162071741e21   4.72236648286965e21
      54   -3.54177486215224e21   1.17260394005318       2.36118324143483e21
      55   -1.77088743107612e21   1.17260394005318       1.18059162071742e21
      56   -8.85443715538059e20   1.17260394005318       5.90295810358706e20
      57   -4.42721857769030e20   1.17260394005318       3.68934881474192e20
      58   -2.21360928884515e20   3.68934881474191e19    1.47573952589677e20
      59   -9.22337203685478e19   -1.84467440737096e19   9.22337203685478e19
      60   -5.53402322211287e19   -1.84467440737096e19   4.61168601842739e19
      61   -3.68934881474192e19   4.61168601842739e18    1.38350580552822e19
      62   -1.38350580552822e19   -2.30584300921369e18   9.22337203685478e18
      63   -8.07045053224793e18   1.15292150460685e18    3.45876451382055e18
      64   -3.45876451382055e18   5.76460752303423e17    2.30584300921370e18
      65   -1.44115188075856e18   1.17260394005318       1.44115188075856e18
      66   -7.20575940379280e17   1.44115188075856e17    7.20575940379280e17
      67   -2.88230376151712e17   1.17260394005318       3.60287970189640e17
      68   -2.16172782113784e17   -7.20575940379279e16   1.80143985094820e17
      69   -1.26100789566374e17   1.80143985094820e16    5.40431955284460e16
      70   -6.30503947831870e16   1.17260394005318       3.60287970189640e16
      71   -2.70215977642230e16   -4.50359962737049e15   1.80143985094820e16
      72   -1.57625986957968e16   2.25179981368525e15    6.75539944105575e15
      73   -6.75539944105575e15   1.12589990684263e15    2.25179981368525e15
      74   -2.81474976710656e15   1.17260394005318       1.68884986026394e15
      75   -1.40737488355328e15   -2.81474976710655e14   8.44424930131970e14
      76   -8.44424930131967e14   1.17260394005318       4.22212465065986e14
      77   -4.22212465065983e14   1.17260394005318       2.11106232532994e14
      78   -2.11106232532991e14   1.17260394005318       7.03687441776652e13
      79   -7.03687441776629e13   1.75921860444172e13    3.51843720888332e13
      80   -4.39804651110389e13   1.17260394005318       1.75921860444172e13
      81   -1.75921860444149e13   1.17260394005318       8.79609302220918e12
      82   -1.09951162777589e13   1.17260394005318       4.39804651110518e12
      83   -4.39804651110283e12   1.17260394005318       3.29853488332918e12
      84   -2.19902325555083e12   5.49755813889173e11    2.19902325555318e12
      85   -8.24633720830828e11   2.74877906945173e11    8.24633720833173e11
      86   -2.74877906942828e11   1.17260394005318       5.49755813889173e11
      87   -1.37438953470828e11   1.17260394005318       2.74877906945173e11
      88   -6.87194767348274e10   1.17260394005318       1.37438953473173e11
      89   -3.43597383668274e10   1.17260394005318       6.87194767371727e10
      90   -1.71798691828274e10   8.58993459317260e9     1.71798691851727e10
      91   -4.29496729482740e9    -4.29496729482740e9    1.28849018891727e10
      92   -4.29496729482740e9    1.17260394005318       6.44245094517261e9
      93   -2.14748364682740e9    1.07374182517260e9     2.14748364917261e9
      94   -5.36870910827397e8    1.17260394005318       1.61061273717261e9
      95   -2.68435454827397e8    -2.68435454827396e8    8.05306369172604e8
      96   -2.68435454827397e8    1.17260394005318       2.68435457172604e8
      97   -6.71088628273961e7    -6.71088628273961e7    1.34217729172604e8
      98   -6.71088628273961e7    3.35544331726039e7     3.35544331726040e7
      99   -1.67772148273961e7    1.17260394005318       3.35544331726040e7
     100   -8.38860682739606e6    1.17260394005318       8.38860917260395e6
     101   -4.19430282739606e6    1.17260394005318       4.19430517260395e6
     102   -2.09715082739606e6    1.17260394005318       2.09715317260395e6
     103   -1.04857482739606e6    1.17260394005318       1.04857717260395e6
     104   -5.24286827396060e5    1.17260394005318       5.24289172603941e5
     105   -2.62142827396060e5    1.17260394005318       2.62145172603941e5
     106   -1.31070827396060e5    1.17260394005318       1.31073172603941e5
     107   -6.55348273960600e4    1.17260394005318       1.17260394005318
     108   -3.27668273960600e4    1.17260394005318       1.17260394005318
     109   -1.63828273960600e4    1.17260394005318       1.17260394005318
     110   -8.19082739605995e3    1.17260394005318       1.17260394005318
     111   -4.09482739605995e3    1.17260394005318       1.17260394005318
     112   -2.04682739605995e3    1.17260394005318       1.17260394005318
     113   -1.02282739605995e3    1.17260394005318       1.17260394005318
     114   -5.10827396059947e2    1.17260394005318       1.17260394005318
     115   -2.54827396059947e2    1.17260394005318       1.17260394005318
     116   -1.26827396059947e2    1.17260394005318       1.17260394005318
     117   -6.28273960599469e1    1.17260394005318       1.17260394005318
     118   -3.08273960599469e1    1.17260394005318       1.17260394005318
     119   -1.48273960599469e1    1.17260394005318       1.17260394005318
     120   -6.82739605994683      1.17260394005318       1.17260394005318
     121   -2.82739605994683      1.17260394005318       1.17260394005318
     122   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
     123   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
     124   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
     125   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
     126   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
     127   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
     128   -8.27396059946822e-1   -8.27396059946821e-1   -8.27396059946821e-1
*/