/* $Id: pb-poly.c 47490 2011-11-10 17:07:08Z vinc17/ypig $
 *
 * Written during CNC'2 (June 2009). Needs MPFR 2.4 for mpfr_printf.
 */

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

static void poly (mpfr_t y, mpfr_t x, int d, mpfr_t *a, int snflag)
{
  int inex;

  mpfr_set_ui (y, 0, MPFR_RNDN);
  while (d >= 0)
    {
      inex = mpfr_mul (y, y, x, MPFR_RNDN);
      if (snflag)
        mpfr_subnormalize (y, inex, MPFR_RNDN);
      inex = mpfr_add (y, y, a[d], MPFR_RNDN);
      if (snflag)
        mpfr_subnormalize (y, inex, MPFR_RNDN);
      d--;
    }
}

int main (int argc, char *argv[])
{
  mpfr_t *mp, my;
  double *fp, fy;
  int d, i;

  if (argc < 3)
    {
      fprintf (stderr, "Usage: pb-poly <x> <a0> ...\n");
      exit (EXIT_FAILURE);
    }
  d = argc - 3;

  mp = malloc ((d + 2) * sizeof (*mp));
  if (mp == NULL)
    {
      fprintf (stderr, "pb-poly: malloc failed (mpfr_t)\n");
      exit (EXIT_FAILURE);
    }

  fp = malloc ((d + 2) * sizeof (*fp));
  if (fp == NULL)
    {
      fprintf (stderr, "pb-poly: malloc failed (double)\n");
      exit (EXIT_FAILURE);
    }

  for (i = 0; i <= d + 1; i++)
    {
      mpfr_init2 (mp[i], 53);
#if 1
      fp[i] = atof (argv[i+1]);
      mpfr_set_d (mp[i], fp[i], MPFR_RNDN);
#else
      if (mpfr_set_str (mp[i], argv[i+1], 0, MPFR_RNDN))
        {
          fprintf (stderr, "pb-poly: bad argument %i (%s)\n",
                   i+1, argv[i+1]);
          exit (EXIT_FAILURE);
        }
      fp[i] = mpfr_get_d (mp[i], MPFR_RNDN);
#endif
    }

  mpfr_init2 (my, 53);
  poly (my, mp[0], d, mp + 1, 0);
  mpfr_printf ("my = %.17Rg\n", my);

  mpfr_set_emin (-1073);
  mpfr_set_emax (1024);
  for (i = 0; i <= d + 1; i++)
    if (mpfr_number_p (mp[i]) && ! mpfr_zero_p (mp[i]))
      {
        mp_exp_t e;
        e = mpfr_get_exp (mp[i]);
        assert (e >= -1073 && e <= 1024);
      }
  poly (my, mp[0], d, mp + 1, 1);
  mpfr_printf ("my = %.17Rg (with mpfr_subnormalize)\n", my);

  fy = 0;
  for (i = d + 1; i > 0; i--)
    fy = fy * fp[0] + fp[i];
  printf ("fy = %.17g\n", fy);

  return 0;
}