xref: /dports/math/mpsolve/mpsolve-3.2.1/src/libmps/monomial/monomial-poly.c

/*
 * This file is part of MPSolve 3.2.1
 *
 * Copyright (C) 2001-2020, Dipartimento di Matematica "L. Tonelli", Pisa.
 * License: http://www.gnu.org/licenses/gpl.html GPL version 3 or higher
 *
 * Authors:
 *   Dario Andrea Bini <bini@dm.unipi.it>
 *   Giuseppe Fiorentino <fiorent@dm.unipi.it>
 *   Leonardo Robol <leonardo.robol@unipi.it>
 */

#include <mps/mps.h>
#include <string.h>
#include <assert.h>

/**
 * @brief Return a newly allocated mps_monomial_poly of the given degree.
 */
mps_monomial_poly *
mps_monomial_poly_new (mps_context * s, long int degree)
{
  int i;
  mps_monomial_poly  * mp = mps_new (mps_monomial_poly);

  mps_polynomial_init (s, MPS_POLYNOMIAL (mp));

  /* Load monomial-poly methods */
  mps_polynomial *poly = (mps_polynomial*)mp;
  poly->type_name = "mps_monomial_poly";
  poly->feval = mps_monomial_poly_feval;
  poly->deval = mps_monomial_poly_deval;
  poly->meval = mps_monomial_poly_meval;
  poly->fstart = mps_monomial_poly_fstart;
  poly->dstart = mps_monomial_poly_dstart;
  poly->mstart = mps_monomial_poly_mstart;
  poly->free = mps_monomial_poly_free;
  poly->raise_data = mps_monomial_poly_raise_precision;
  poly->fnewton = mps_monomial_poly_fnewton;
  poly->dnewton = mps_monomial_poly_dnewton;
  poly->mnewton = mps_monomial_poly_mnewton;
  poly->get_leading_coefficient = mps_monomial_poly_get_leading_coefficient;

  /* Set the degree of the polynomial */
  MPS_POLYNOMIAL (mp)->degree = degree;

  /* Allocate the space for the coefficients of the polynomial, all
   * the floating point versions. */
  mp->spar = mps_boolean_valloc (degree + 2);
  mp->fpc = cplx_valloc (degree + 1);
  mp->fpr = double_valloc (degree + 1);
  mp->dpr = rdpe_valloc (degree + 1);
  mp->dpc = cdpe_valloc (degree + 1);
  mp->db.mfpc1 = mpc_valloc (degree + 1);
  mp->db.mfpc2 = mpc_valloc (degree + 1);
  mp->mfpc = mp->db.mfpc1;
  mp->db.active = 1;
  mp->mfpr = mpf_valloc (degree + 1);

  mpf_vinit2 (mp->mfpr, degree + 1, s->mpwp);
  mpc_vinit2 (mp->db.mfpc1, degree + 1, s->mpwp);
  mpc_vinit2 (mp->db.mfpc2, degree + 1, s->mpwp);

  /* Allocate space for the moduli of the coefficients */
  mp->fap = double_valloc (degree + 1);
  mp->dap = rdpe_valloc (degree + 1);

  /* Allocate space for the coefficients of the derivative */
  mp->fppc = cplx_valloc (degree + 1);
  mp->mfppc = mpc_valloc (degree + 1);
  mpc_vinit2 (mp->mfppc, degree + 1, s->mpwp);

  /* Allocate space for the coefficients initially parsed as
   * exact */
  mp->initial_mqp_r = mpq_valloc (degree + 1);
  mp->initial_mqp_i = mpq_valloc (degree + 1);

  mpq_vinit (mp->initial_mqp_r, degree + 1);
  mpq_vinit (mp->initial_mqp_i, degree + 1);

  pthread_mutex_init (&mp->regenerating, NULL);

  mp->mfpc_mutex = mps_newv (pthread_mutex_t, degree + 1);
  for (i = 0; i <= degree; i++)
    pthread_mutex_init (&mp->mfpc_mutex[i], NULL);

  memset (mp->spar, 0, sizeof(mps_boolean) * (degree + 1));
  memset (mp->fap, 0, sizeof(double) * (degree + 1));

  for (i = 0; i <= degree; i++)
    rdpe_set (mp->dap[i], rdpe_zero);

  MPS_POLYNOMIAL (mp)->structure = MPS_STRUCTURE_UNKNOWN;
  mp->prec = s->mpwp;

  return mp;
}

/**
 * @brief Free a instance of <code>mps_monomial_poly</code> previously
 * allocated with <code>mps_monomial_poly_new()</code>.
 */
void
mps_monomial_poly_free (mps_context * s, mps_polynomial * p)
{
  mps_monomial_poly *mp = MPS_MONOMIAL_POLY (p);

  mps_boolean_vfree (mp->spar);
  double_vfree (mp->fpr);
  cplx_vfree (mp->fpc);
  rdpe_vfree (mp->dpr);
  cdpe_vfree (mp->dpc);

  double_vfree (mp->fap);
  rdpe_vfree (mp->dap);

  mpf_vclear (mp->mfpr, MPS_POLYNOMIAL (mp)->degree + 1);
  mpc_vclear (mp->db.mfpc1, MPS_POLYNOMIAL (mp)->degree + 1);
  mpc_vclear (mp->db.mfpc2, MPS_POLYNOMIAL (mp)->degree + 1);

  mpf_vfree (mp->mfpr);
  mpc_vfree (mp->db.mfpc1);
  mpc_vfree (mp->db.mfpc2);

  mpq_vclear (mp->initial_mqp_r, MPS_POLYNOMIAL (mp)->degree + 1);
  mpq_vclear (mp->initial_mqp_i, MPS_POLYNOMIAL (mp)->degree + 1);

  mpq_vfree (mp->initial_mqp_r);
  mpq_vfree (mp->initial_mqp_i);

  cplx_vfree (mp->fppc);
  mpc_vclear (mp->mfppc, MPS_POLYNOMIAL (mp)->degree + 1);
  mpc_vfree (mp->mfppc);

  free (mp->mfpc_mutex);

  free (mp);
}

/**
 * @brief Raise the precision bits of the multiprecision fields of the
 * polynomial to selected value.
 *
 * If the polynomial coefficients were generated by integer or rational
 * input they will be autoregenerated.
 *
 * @param s  The status of the computation.
 * @param p The polynomial that need the precision raised, casted to a mps_polynomial.
 * @param prec The selected bits of precision.
 * @return The precision set.
 */
long int
mps_monomial_poly_raise_precision (mps_context * s, mps_polynomial * p, long int prec)
{
  mps_monomial_poly *mp = MPS_MONOMIAL_POLY (p);
  int k;
  mpc_t * raising_mfpc;

  pthread_mutex_lock (&mp->regenerating);

  /* For floating point coefficients, it might happen that
   * the coefficients already have higher precision. In that
   * case, we return now. */
  if (prec <= mp->prec ||
      (MPS_STRUCTURE_IS_FP(p->structure) && mpc_get_prec(mp->mfpc[0]) >= prec))
    {
      MPS_DEBUG_WITH_INFO(s, "Not increasing precision, the coefficients are already at the required accuracy");
      pthread_mutex_unlock (&mp->regenerating);
      return mp->prec;
    }

  if (mp->db.active == 1)
    raising_mfpc = mp->db.mfpc2;
  else
    raising_mfpc = mp->db.mfpc1;

  if (MPS_STRUCTURE_IS_FP(p->structure))
    {
      long int current_wp = mpc_get_prec(mp->mfpc[0]);
      long int raising_wp = mpc_get_prec(raising_mfpc[0]);

      if (raising_wp > current_wp)
        {
          for (k = 0; k <= s->n; k++)
            {
              mpc_set_prec (mp->mfpc[k], raising_wp);
              mpc_set (mp->mfpc[k], raising_mfpc[k]);
            }
        }
    }

  /* raise the precision of  mfpc */
  if (MPS_IS_MONOMIAL_POLY (p))
    for (k = 0; k < MPS_POLYNOMIAL (mp)->degree + 1; k++)
      {
        mpc_set_prec (raising_mfpc[k], prec);
      }

  /* Raise the precision of p' */
  if (MPS_DENSITY_IS_SPARSE (p->density))
    for (k = 0; k < MPS_POLYNOMIAL (mp)->degree; k++)
      if (mp->spar[k + 1])
        {
          mpc_set_prec (mp->mfppc[k], prec);
          mpc_mul_ui (mp->mfppc[k], mp->mfpc[k + 1], k + 1);
        }

  if (MPS_STRUCTURE_IS_INTEGER (p->structure)
      || MPS_STRUCTURE_IS_RATIONAL (p->structure))
    {
      for (k = 0; k <= MPS_POLYNOMIAL (mp)->degree; ++k)
        {
          mpf_set_q (mpc_Re (raising_mfpc[k]), mp->initial_mqp_r[k]);
          mpf_set_q (mpc_Im (raising_mfpc[k]), mp->initial_mqp_i[k]);
          mpc_rmod (mp->dap[k], raising_mfpc[k]);
          mp->fap[k] = rdpe_get_d (mp->dap[k]);
        }
    }
  else
    {
      for (k = 0; k <= MPS_POLYNOMIAL (mp)->degree; k++)
        {
          mpc_set (raising_mfpc[k], mp->mfpc[k]);
        }
    }

  mp->db.active = (mp->db.active % 2) + 1;
  mp->mfpc = raising_mfpc;

  mp->prec = prec;
  pthread_mutex_unlock (&mp->regenerating);

  return mp->prec;
}

long int
mps_monomial_poly_get_precision (mps_context * s, mps_monomial_poly * mp)
{
  return mpc_get_prec (mp->mfpc[0]);
}

/**
 * @brief This routine can be used to set the i-th coefficients of the
 * polynomial with a multiprecision rational number.
 *
 * @param s The mps_context the will be used for the computation. This shall be passed
 * along with the polynomial because it manages user interaction and can perform
 * error handling.
 *
 * @param mp The monomial_poly that will hold the coefficient.
 * @param i  The exponent of the monomial to which the coefficient is related.
 * @param real_part The real part of the coefficient.
 * @param imag_part The imaginary part of the coefficients.
 */
void
mps_monomial_poly_set_coefficient_q (mps_context * s, mps_monomial_poly * mp, long int i,
                                     mpq_t real_part, mpq_t imag_part)
{
  /* Updating data_type information */
  if (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_UNKNOWN)
    MPS_POLYNOMIAL (mp)->structure = (mpq_sgn (imag_part) != 0) ?
                                     MPS_STRUCTURE_COMPLEX_RATIONAL : MPS_STRUCTURE_REAL_RATIONAL;

  if (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_RATIONAL &&
      mpq_sgn (imag_part) != 0)
    MPS_POLYNOMIAL (mp)->structure = MPS_STRUCTURE_COMPLEX_RATIONAL;

  assert (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_COMPLEX_RATIONAL ||
          MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_RATIONAL ||
          MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_COMPLEX_INTEGER ||
          MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_INTEGER);

  /* Set the coefficients first */
  mpq_set (mp->initial_mqp_r[i], real_part);
  mpq_set (mp->initial_mqp_i[i], imag_part);

  /* Then update the other coefficients */
  mpf_set_q (mpc_Re (mp->mfpc[i]), real_part);
  mpf_set_q (mpc_Im (mp->mfpc[i]), imag_part);

  mpc_get_cdpe (mp->dpc[i], mp->mfpc[i]);
  mpc_get_cplx (mp->fpc[i], mp->mfpc[i]);

  if ((mpq_sgn (real_part) == 0) && (mpq_sgn (imag_part) == 0))
    mp->spar[i] = false;
  else
    mp->spar[i] = true;

  if (mp->spar[i])
    {
      mpc_get_cplx (mp->fpc[i], mp->mfpc[i]);
      mpc_get_cdpe (mp->dpc[i], mp->mfpc[i]);

      /* Compute modules of coefficients */
      cdpe_mod (mp->dap[i], mp->dpc[i]);
      mp->fap[i] = rdpe_get_d (mp->dap[i]);

      if (i > 0)
        mpc_mul_ui (mp->mfppc[i - 1], mp->mfppc[i], i);
    }
  else
    {
      cplx_set (mp->fpc[i], cplx_zero);
      cdpe_set (mp->dpc[i], cdpe_zero);

      rdpe_set (mp->dap[i], rdpe_zero);
      mp->fap[i] = 0.0f;
    }
}

/**
 * @brief Set the coefficient in position i of the mpnomial.
 *
 * @param s The <code>mps_context</code> associated to this computation.
 * @param mp The <code>monomial_poly</code> in which the coefficients will be set.
 * @param i The index of the coefficient to set.
 * @param real_part The real part of the coefficient.
 * @param imag_part The imaginary part of the coefficient.
 */
void
mps_monomial_poly_set_coefficient_d (mps_context * s, mps_monomial_poly * mp, long int i,
                                     double real_part, double imag_part)
{
  /* Updating data structure information */
  if (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_UNKNOWN)
    MPS_POLYNOMIAL (mp)->structure = (imag_part == 0) ?
                                     MPS_STRUCTURE_REAL_FP : MPS_STRUCTURE_COMPLEX_FP;

  if (imag_part != 0 && MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_FP)
    MPS_POLYNOMIAL (mp)->structure = MPS_STRUCTURE_COMPLEX_FP;

  assert (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_FP ||
          MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_COMPLEX_FP);

  /* Set the coefficients */
  mpf_set_d (mp->mfpr[i], real_part);
  mpc_set_d (mp->mfpc[i], real_part, imag_part);

  /* Update spar */
  mp->spar[i] = !((real_part == 0) && (imag_part == 0));

  if (mp->spar[i])
    {
      mpc_get_cplx (mp->fpc[i], mp->mfpc[i]);
      mpc_get_cdpe (mp->dpc[i], mp->mfpc[i]);

      /* Compute modules of coefficients */
      cdpe_mod (mp->dap[i], mp->dpc[i]);
      mp->fap[i] = rdpe_get_d (mp->dap[i]);

      if (i > 0)
        mpc_mul_ui (mp->mfppc[i - 1], mp->mfppc[i], i);
    }
  else
    {
      cplx_set (mp->fpc[i], cplx_zero);
      cdpe_set (mp->dpc[i], cdpe_zero);

      rdpe_set (mp->dap[i], rdpe_zero);
      mp->fap[i] = 0.0f;
    }
}

void
mps_monomial_poly_set_coefficient_int (mps_context * s, mps_monomial_poly * mp, long int i,
                                       long long real_part, long long imag_part)
{
  /* Updating data_type information */
  if (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_UNKNOWN)
    MPS_POLYNOMIAL (mp)->structure = (imag_part) != 0 ?
                                     MPS_STRUCTURE_COMPLEX_INTEGER : MPS_STRUCTURE_REAL_INTEGER;

  if (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_INTEGER &&
      imag_part != 0)
    MPS_POLYNOMIAL (mp)->structure = MPS_STRUCTURE_COMPLEX_INTEGER;

  assert (MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_COMPLEX_INTEGER ||
          MPS_POLYNOMIAL (mp)->structure == MPS_STRUCTURE_REAL_INTEGER);

  /* Set the coefficients first */
  mpq_set_si (mp->initial_mqp_r[i], real_part, 1);
  mpq_set_si (mp->initial_mqp_i[i], imag_part, 1);

  /* Then update the other coefficients */
  mpf_set_q (mpc_Re (mp->mfpc[i]), mp->initial_mqp_r[i]);
  mpf_set_q (mpc_Im (mp->mfpc[i]), mp->initial_mqp_i[i]);

  mpc_get_cdpe (mp->dpc[i], mp->mfpc[i]);
  mpc_get_cplx (mp->fpc[i], mp->mfpc[i]);

  if ((real_part == 0) && (imag_part == 0))
    mp->spar[i] = false;
  else
    mp->spar[i] = true;

  if (mp->spar[i])
    {
      mpc_get_cplx (mp->fpc[i], mp->mfpc[i]);
      mpc_get_cdpe (mp->dpc[i], mp->mfpc[i]);

      /* Compute modules of coefficients */
      cdpe_mod (mp->dap[i], mp->dpc[i]);
      mp->fap[i] = rdpe_get_d (mp->dap[i]);

      if (i > 0)
        mpc_mul_ui (mp->mfppc[i - 1], mp->mfppc[i], i);
    }
  else
    {
      cplx_set (mp->fpc[i], cplx_zero);
      cdpe_set (mp->dpc[i], cdpe_zero);

      rdpe_set (mp->dap[i], rdpe_zero);
      mp->fap[i] = 0.0f;
    }
}

void
mps_monomial_poly_set_coefficient_f (mps_context * s, mps_monomial_poly * p, long int i,
                                     mpc_t coeff)
{
  /* Updating data_type information */
  if (MPS_POLYNOMIAL (p)->structure == MPS_STRUCTURE_UNKNOWN)
    MPS_POLYNOMIAL (p)->structure = MPS_STRUCTURE_COMPLEX_FP;

  mpc_set_prec (p->mfpc[i], mpc_get_prec(coeff));
  mpc_set (p->mfpc[i], coeff);
  mpc_get_cplx (p->fpc[i], p->mfpc[i]);
  mpc_get_cdpe (p->dpc[i], p->mfpc[i]);

  mpc_rmod (p->dap[i], p->mfpc[i]);
  p->fap[i] = rdpe_get_d (p->dap[i]);

  p->spar[i] = !mpc_eq_zero (coeff);

  if (i > 0)
    mpc_mul_ui (p->mfppc[i - 1], p->mfppc[i], i);
}

/**
 * @brief Set the \f$i\f$-th coefficient of the polynomial.
 *
 * @param s The current mps_context.
 * @param p The polynomial whose coefficient should be set.
 * @param i The degree of the coefficient to set.
 * @param real_coeff A string representing the real part of the coefficient
 * that should be set, or NULL as an alias for \f$0\f$.
 * @param imag_coeff A string representing the imaginary part of the coefficient
 * that should be set, or NULL as an alias for \f$0\f$.
 */
void
mps_monomial_poly_set_coefficient_s (mps_context * s, mps_monomial_poly * p,
                                     int i, const char * real_coeff,
                                     const char * imag_coeff)
{
  char * eq_real = mps_utils_build_equivalent_rational_string (s, real_coeff);
  char * eq_imag = mps_utils_build_equivalent_rational_string (s, imag_coeff);

  mpq_t tmp_real, tmp_imag;

  mpq_init (tmp_real);
  mpq_init (tmp_imag);

  if (eq_real != NULL)
    mpq_set_str (tmp_real, eq_real, 10);
  else
    mps_warn (s, "Invalid input for mps_monomial_set_coefficient_s: %s", real_coeff);

  if (eq_imag != NULL)
    mpq_set_str (tmp_imag, eq_imag, 10);
  else
    mps_warn (s, "Invalid input for mps_monomial_set_coefficient_s: %s", imag_coeff);

  mps_monomial_poly_set_coefficient_q (s, p, i, tmp_real, tmp_imag);

  mpq_clear (tmp_real);
  mpq_clear (tmp_imag);

  free (eq_real);
  free (eq_imag);
}

/**
 * @brief Get a double version of the \f$i\f$-th coefficient of the polynomial.
 *
 * @param s The current mps_context.
 * @param p The polynomial of which you want to know the coefficient.
 * @param i The degree of the coefficient to obtain.
 * @param output A cplx_t where the result will be stored.
 */
void mps_monomial_poly_get_coefficient_d (mps_context * s, mps_monomial_poly * p,
                                          int i, cplx_t output)
{
  if (i >= 0 && i <= MPS_POLYNOMIAL (p)->degree)
    cplx_set (output, p->fpc[i]);
  else
    cplx_set (output, cplx_zero);
}

/**
 * @brief Get a rational version of the \f$i\f$-th coefficient of the polynomial.
 *
 * @param s The current mps_context.
 * @param p The polynomial of which you want to know the coefficient.
 * @param i The degree of the coefficient to obtain.
 * @param real_output A mpq_t where the real part of the result will be stored.
 * @param imag_output A mpq_t where the imaginary part of the result will be stored.
 */
void mps_monomial_poly_get_coefficient_q (mps_context * s, mps_monomial_poly * p,
                                          int i, mpq_t real_output, mpq_t imag_output)
{
  mps_polynomial *poly = MPS_POLYNOMIAL (p);

  if ((!MPS_STRUCTURE_IS_RATIONAL (poly->structure)) &&
      (!MPS_STRUCTURE_IS_INTEGER (poly->structure)))
    {
      mps_error (s, "Cannot extract rational coefficients from a floating point polynomial");
      return;
    }

  if (i >= 0 && i <= poly->degree)
    {
      mpq_set (real_output, p->initial_mqp_r[i]);
      mpq_set (imag_output, p->initial_mqp_i[i]);
    }
  else
    {
      mpq_set_si (real_output, 0, 1);
      mpq_set_si (imag_output, 0, 1);
    }
}

/**
 * @brief Get the k-th derivative of p with floating point coefficients
 * approximated with the precision wp.
 *
 * @param s The current mps_context
 * @param p The polynomial to derive
 * @param k The order of the derivative that should be computed.
 * @param wp The selected output working precision
 */
mps_monomial_poly *
mps_monomial_poly_derive (mps_context * s, mps_monomial_poly * p, int k, long int wp)
{
  mps_monomial_poly * d = mps_monomial_poly_new (s, MPS_POLYNOMIAL (p)->degree - k);
  int i, j;

  MPS_POLYNOMIAL (d)->structure = MPS_POLYNOMIAL (p)->structure;
  MPS_POLYNOMIAL (d)->density = MPS_POLYNOMIAL (p)->density;
  MPS_POLYNOMIAL (d)->prec = MPS_POLYNOMIAL (p)->prec;

  if (wp != s->mpwp)
    mps_monomial_poly_raise_precision (s, MPS_POLYNOMIAL (d), wp);

  switch (MPS_POLYNOMIAL (p)->structure)
    {
    case MPS_STRUCTURE_COMPLEX_INTEGER:
    case MPS_STRUCTURE_COMPLEX_RATIONAL:
    case MPS_STRUCTURE_REAL_INTEGER:
    case MPS_STRUCTURE_REAL_RATIONAL:
    {
      mpq_t coeff_r;
      mpq_t coeff_i;
      mpq_t qtmp;
      mpq_init (coeff_r);
      mpq_init (coeff_i);
      mpq_init (qtmp);

      for (i = 0; i <= MPS_POLYNOMIAL (d)->degree; i++)
        {
          mpq_set (coeff_r, p->initial_mqp_r[i + k]);
          mpq_set (coeff_i, p->initial_mqp_i[i + k]);

          for (j = 0; j < k; j++)
            {
              mpq_set_si (qtmp, k + i - j, 1);
              mpq_mul (coeff_r, coeff_r, qtmp);
              mpq_mul (coeff_i, coeff_i, qtmp);
            }
          mps_monomial_poly_set_coefficient_q (s, d, i, coeff_r, coeff_i);
        }

      mpq_clear (coeff_r);
      mpq_clear (coeff_i);
      mpq_clear (qtmp);
    }

    break;

    default:
    {
      mpc_t * coeffs = mps_newv (mpc_t, MPS_POLYNOMIAL (d)->degree + 1);
      mpc_vinit2 (coeffs, MPS_POLYNOMIAL (d)->degree + 1, wp);

      for (i = 0; i <= MPS_POLYNOMIAL (d)->degree; i++)
        mpc_set (coeffs[i], p->mfpc[i + k]);

      for (j = 0; j < k; j++)
        for (i = 0; i <= MPS_POLYNOMIAL (d)->degree; i++)
          mpc_mul_eq_ui (coeffs[i], k + i - j);


      for (i = 0; i <= MPS_POLYNOMIAL (d)->degree; i++)
        mps_monomial_poly_set_coefficient_f (s, d, i, coeffs[i]);

      mpc_vclear (coeffs, MPS_POLYNOMIAL (d)->degree + 1);
      free (coeffs);
    }
    break;
    }

  if (MPS_DENSITY_IS_SPARSE (MPS_POLYNOMIAL (d)->density))
    for (i = 0; i < MPS_POLYNOMIAL (d)->degree; i++)
      mpc_mul_ui (d->mfppc[i], d->mfpc[i + 1], i + 1);

  return d;
}

mps_boolean
mps_monomial_poly_feval (mps_context *ctx, mps_polynomial *p, cplx_t x, cplx_t value, double * error)
{
  mps_fhorner_with_error (ctx, (mps_monomial_poly*)p, x, value, error);
  return true;
}

mps_boolean
mps_monomial_poly_deval (mps_context *ctx, mps_polynomial *p, cdpe_t x, cdpe_t value, rdpe_t error)
{
  mps_dhorner_with_error (ctx, (mps_monomial_poly*)p, x, value, error);
  return true;
}

mps_boolean
mps_monomial_poly_meval (mps_context * ctx, mps_polynomial *p, mpc_t x, mpc_t value, rdpe_t error)
{
  mps_mhorner_with_error2 (ctx, (mps_monomial_poly*)p, x, value, error, mpc_get_prec (x));
  return true;
}

void
mps_monomial_poly_fstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations)
{
  mps_monomial_poly *mp = MPS_MONOMIAL_POLY (p);

  mps_fstart (ctx, ctx->n, NULL, 0.0, 0.0, ctx->eps_out, mp->fap);
}

void
mps_monomial_poly_dstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations)
{
  mps_monomial_poly *mp = MPS_MONOMIAL_POLY (p);

  mps_dstart (ctx, ctx->n, NULL, (__rdpe_struct*)rdpe_zero,
              (__rdpe_struct*)rdpe_zero, ctx->eps_out,
              mp->dap);
}

void
mps_monomial_poly_mstart (mps_context * ctx, mps_polynomial * p, mps_approximation ** approximations)
{
  mps_monomial_poly *mp = MPS_MONOMIAL_POLY (p);
  mpc_t gg;

  mpc_init2 (gg, ctx->mpwp);
  mpc_set_ui (gg, 0U, 0U);
  mps_mstart (ctx, ctx->n, NULL, (__rdpe_struct*)rdpe_zero,
              (__rdpe_struct*)rdpe_zero, mp->dap, gg);
  mpc_clear (gg);
}

void
mps_monomial_poly_fnewton (mps_context * ctx, mps_polynomial * p,
                           mps_approximation * root, cplx_t corr)
{
  mps_fnewton (ctx, p, root, corr);
}

void
mps_monomial_poly_dnewton (mps_context * ctx, mps_polynomial * p,
                           mps_approximation * root, cdpe_t corr)
{
  mps_dnewton (ctx, p, root, corr);
}

void
mps_monomial_poly_mnewton (mps_context * ctx, mps_polynomial * p,
                           mps_approximation * root, mpc_t corr, long int wp)
{
  mps_mnewton (ctx, p, root, corr, wp);
}

void
mps_monomial_poly_get_leading_coefficient (mps_context * ctx, mps_polynomial * p,
                                           mpc_t leading_coefficient)
{
  mpc_set (leading_coefficient, MPS_MONOMIAL_POLY (p)->mfpc[p->degree]);
}

void mps_monomial_poly_deflate (mps_context * ctx, mps_polynomial * poly)
{
  int i;
  mps_monomial_poly * p = MPS_MONOMIAL_POLY (poly);

  /* count number of zero roots (undeflated input polynomial) */
  int zero_roots = 0;

  while (rdpe_eq (p->dap[zero_roots], rdpe_zero) && zero_roots < poly->degree)
    zero_roots++;

  /* shift down input vectors */
  if (zero_roots)
    {
      for (i = 0; i <= poly->degree - zero_roots; i++)
        {
          rdpe_set (p->dap[i], p->dap[i + zero_roots]);
          p->fap[i] = p->fap[i + zero_roots];
          p->fpr[i] = p->fpr[i + zero_roots];
          cplx_set (p->fpc[i], p->fpc[i + zero_roots]);
          rdpe_set (p->dpr[i], p->dpr[i + zero_roots]);
          cdpe_set (p->dpc[i], p->dpc[i + zero_roots]);
          mpf_set (p->mfpr[i], p->mfpr[i + zero_roots]);
          mpc_set (p->mfpc[i], p->mfpc[i + zero_roots]);
          if (i < poly->degree - zero_roots)
            mpc_set (p->mfppc[i], p->mfppc[i + zero_roots]);
          mpq_set (p->initial_mqp_r[i], p->initial_mqp_r[i + zero_roots]);
          mpq_set (p->initial_mqp_i[i], p->initial_mqp_i[i + zero_roots]);
          p->spar[i] = p->spar[i + zero_roots];
        }
    }

  /* FIXME: We need to reallocate the correct storage for the
   * polynomial. */

  poly->degree -= zero_roots;
}