/* Ecole Calcul Numerique Certifie 2
   Nancy - 25-26 juin 2009
   probleme 1.12 */

#include <stdlib.h>
#include <stdio.h>
#include <mpc.h>

typedef struct
{
  mp_prec_t prec;
  int degree;
  mpc_t *coeff;
} poly_t;

void
polynomial_init (poly_t *p, int degree, mp_prec_t prec)
{
  int d;

  p->prec = prec;
  p->degree = degree;
  p->coeff = (mpc_t *) malloc ((degree+1)*sizeof(mpc_t));
  for (d = 0; d <= degree; d++)
    {
      mpc_init2 (p->coeff[d], prec);
      mpc_set_ui_ui (p->coeff[d], 0, 0, MPC_RNDNN);
    }
}

void
polynomial_clear (poly_t *p)
{
  int d;
  for (d = 0; d <= p->degree; d++)
    mpc_clear (p->coeff[d]);
  free (p->coeff);
}

int
polynomial_read_coeff (poly_t *p, FILE *stream)
{
  int d;

  for (d = p->degree; d >= 0; d--)
    {
      if (stream == NULL)
        printf ("coefficient [%d] = ", d);

      if (mpc_inp_str (p->coeff[d], stream, NULL, 10, MPC_RNDNN)==-1)
        {
          fprintf (stderr, "Cannot read mpc value\n");
          return 1;
        }
    }
  return 0;
}

void
polynomial_out_str (FILE *stream, int base, int n, poly_t *p, mpc_rnd_t rnd)
{
  int d;

  if (stream == NULL)
    stream = stdout;

  for (d = p->degree; d >= 0; d--)
    {
      mpc_out_str (stream, base, n, p->coeff[d], rnd);
      if (d > 0)
        fprintf (stream, " * z^%d + ", d);
      else
        fprintf (stream, "\n");
    }
}

void
polynomial_eval (mpc_ptr rop, poly_t *p, mpc_srcptr op, mpc_rnd_t rnd)
{
  mp_prec_t pr, pi;
  mpc_t z;
  int d;

  mpc_get_prec2 (&pr, &pi, rop);
  mpc_init3 (z, pr, pi);

  d = p->degree;
  mpc_set (z, p->coeff[d--], rnd);
  for (; d >= 0; d--)
    {
      mpc_mul (z, z, op, rnd);
      mpc_add (z, z, p->coeff[d], rnd);
    }

  mpc_set (rop, z, rnd);
  mpc_clear (z);
}

int
polynomial_mul (poly_t *rop, poly_t *op1, poly_t *op2, mpc_rnd_t rnd)
{
  int i, j, jmin, jmax;
  poly_t *pmin;
  poly_t *pmax;
  mpc_t z;

  if (rop == op1 || rop == op2)
    {
      printf ("Input operands must differ from return operand in polynomial_mul\n");
      return 1;
    }

  if (rop->degree < op1->degree + op2-> degree)
    {
      printf ("Degrees do not agree\n");
      return 1;
    }

  if (op1->degree < op2->degree)
    {
      pmin = op1;
      pmax = op2;
    }
  else
    {
      pmin = op2;
      pmax = op1;
    }

  mpc_init2 (z, rop->prec);
  for (i = 0; i <= rop->degree; i++)
    {
      mpc_set_ui_ui (rop->coeff[i], 0, 0, MPC_RNDNN);

      jmax = i < pmin->degree ? i : pmin->degree;
      jmin = i < pmax->degree ? 0 : i - pmax->degree;
      for (j = jmin; j <= jmax; j++)
        {
          mpc_mul (z, pmin->coeff[j], pmax->coeff[i-j], rnd);
          mpc_add (rop->coeff[i], rop->coeff[i], z, rnd);
        }
    }

  mpc_clear (z);
  return 0;
}

int
main (int argc, char **argv)
{
  mpc_t z;
  poly_t p1, p2, p3;

  mpc_init2 (z, 53);
  polynomial_init (&p1, 2, 53);
  polynomial_init (&p2, 4, 53);
  polynomial_init (&p3, 6, 53);

  printf ("poly 1 ===\n");
  if (polynomial_read_coeff (&p1, NULL))
    return 1;
  polynomial_out_str (NULL, 10, 3, &p1, MPC_RNDNN);

  printf ("\npoly 2 ===\n");
  if (polynomial_read_coeff (&p2, NULL))
    return 1;
  polynomial_out_str (NULL, 10, 3, &p2, MPC_RNDNN);

  if (polynomial_mul (&p3, &p1, &p2, MPC_RNDNN))
    return 0;
  printf ("\npoly 3 ===\n");
  polynomial_out_str (NULL, 10, 3, &p3, MPC_RNDNN);

  while (1)
    {
      printf ("\nArgument value (give an invalid string to exit): ");
      if (mpc_inp_str (z, NULL, NULL, 10, MPC_RNDNN) == -1)
        break;

      polynomial_eval (z, &p3, z, MPC_RNDNN);
      printf ("Result: ");
      mpc_out_str (NULL, 10, 0, z, MPC_RNDNN);
      putchar ('\n');
    }

  mpc_clear (z);
  polynomial_clear (&p1);
  polynomial_clear (&p2);
  polynomial_clear (&p3);
  return 0;
}