/*
mpcx_rem

computes the remainder of f divided by g

Copyright (C) 2009 Andreas Enge

This file is part of the MPFRCX Library.

The MPFRCX Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The MPFRCX Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the MPFRCX library; see the file COPYING.LESSER.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA.
*/

#include "mpfrcx-impl.h"

static void mpcx_rem_naive (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g);
static void mpcx_rem_newton (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g);
static void mpcx_rev (mpcx_ptr h, mpcx_srcptr f, const int exp);
static void mpcx_smul (mpcx_ptr h, mpcx_ptr f, mpcx_ptr g, const int exp);
static void mpcx_sinv (mpcx_ptr h, mpcx_ptr f, const int exp);

/*****************************************************************************/

void mpcx_rem (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g)

{
   if (f->deg < g->deg)
      mpcx_set (r, f);
   else
      if (g->deg >= 32 && f->deg >= 63)
         mpcx_rem_newton (r, f, g);
      else
         mpcx_rem_naive (r, f, g);
}

/*****************************************************************************/

static void mpcx_rem_naive (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g)
   /* computes the remainder of f divided by g by the naive algorithm for    */
   /* deg f >= deg g                                                         */
{
   mpc_t tmp;
   int i, j;

   if (mpc_cmp_si (g->coeff [g->deg], 1) == 0) {
      mpc_init2 (tmp, mpc_get_prec (f->coeff [0]));
      /* the only case we actually use for multievaluation */
      mpcx_set (r, f);
      for (i = f->deg - g->deg; i >= 0; i--)
         for (j = 0; j < g->deg; j++) {
            mpc_mul (tmp, r->coeff [g->deg + i], g->coeff [j], MPC_RNDNN);
            mpc_sub (r->coeff [i + j], r->coeff [i + j], tmp, MPC_RNDNN);
         }
      r->deg = g->deg - 1;
      /* To save memory, one might want to do a realloc on r here:
      mpcx_realloc (r, r->deg + 1);
      */
      mpc_clear (tmp);
   }
   else {
      printf ("*** Houston, we have a problem!\n");
      printf ("*** Calling mpcx_rem_naive with a non-monic divisor.\n");
      printf ("*** Go back programming!\n");
      exit (1);
   }
}

/*****************************************************************************/

static void mpcx_rem_newton (mpcx_ptr r, mpcx_srcptr f, mpcx_srcptr g) {
   /* computes the remainder of f divided by g by Newton iterations for      */
   /* deg f >= deg g                                                         */

   mpcx_t q, tmp;

   mpcx_init (q, f->deg - g->deg + 1, r->prec);
   mpcx_init (tmp, f->deg + 1, r->prec);

   mpcx_rev (tmp, f, 0);
   mpcx_rev (q, g, 0);
   mpcx_sinv (q, q, f->deg - g->deg + 1);
   mpcx_smul (q, tmp, q, f->deg - g->deg + 1);
   mpcx_rev (q, q, f->deg - g->deg);

   mpcx_mul (tmp, q, g);
   mpcx_sub (r, f, tmp);
   if (r->deg >= g->deg)
      r->deg = g->deg - 1;

   mpcx_clear (q);
   mpcx_clear (tmp);
}

/*****************************************************************************/

static void mpcx_rev (mpcx_ptr h, mpcx_srcptr f, const int exp) {
   int    overlap, exp_local, i;
   mpcx_t h_local;

   if (exp == 0)
      exp_local = f->deg;
   else
      exp_local = exp;

   if (exp_local < f->deg) {
      printf ("*** Computing a reverse polynomial of too low ");
      printf ("order in 'mpcx_rev'.\n");
      exit (1);
   }

   overlap = (h == f);
   if (overlap)
      mpcx_init (h_local, exp_local + 1, h->prec);
   else {
      mpcx_mv (h_local, h);
      if (h_local->size < exp_local + 1)
         mpcx_realloc (h_local, exp_local + 1);
   }

   h_local->deg = exp_local;
   for (i = 0; i <= f->deg; i++)
      mpc_set (h_local->coeff [exp_local - i], f->coeff [i], MPC_RNDNN);
   for (i = exp_local - f->deg - 1; i >= 0; i--)
      mpc_set_ui (h_local->coeff [i], 0, MPC_RNDNN);

   while (h_local->deg >= 0
         && mpc_cmp_si (h_local->coeff [h_local->deg], 0) == 0)
      h_local->deg--;
   if (overlap)
      mpcx_clear (h);
   mpcx_mv (h, h_local);
}

/*****************************************************************************/

static void mpcx_smul (mpcx_ptr h, mpcx_ptr f, mpcx_ptr g, const int exp) {
   /* computes the short product h = f*g mod X^exp; for the time being,   */
   /* by a complete product and a truncation                              */

   int f_deg, g_deg, h_deg;

   f_deg = f->deg;
   g_deg = g->deg;
   if (f_deg >= exp)
      f->deg = exp - 1;
   if (g_deg >= exp)
      g->deg = exp - 1;
   mpcx_mul (h, f, g);
   h_deg = h->deg;
   f->deg = f_deg;
   g->deg = g_deg;
   h->deg = h_deg;
      /* The previous line is needed if h is one of f or g. */
   if (h->deg >= exp) {
      h->deg = exp - 1;
      while (h->deg >= 0 && mpc_cmp_si (h->coeff [h->deg], 0) == 0)
         h->deg--;
   }
}

/*****************************************************************************/

static void mpcx_sinv (mpcx_ptr h, mpcx_ptr f, const int exp) {
   /* computes h = f^{-1} mod X^exp by Newton iteration */

   int    overlap;
   mpcx_t h_local, tmp;
   int    exp_local;

   overlap = (h == f);
   if (overlap)
      mpcx_init (h_local, exp, h->prec);
   else
      mpcx_mv (h_local, h);
   mpcx_init (tmp, exp, h->prec);

   exp_local = 1;
   h_local->deg = 0;
   mpc_ui_div (h_local->coeff [0], 1, f->coeff [0], MPC_RNDNN);
   while (exp_local < exp) {
      exp_local = (2 * exp_local < exp ? 2 * exp_local : exp);
      mpcx_smul (tmp, h_local, f, exp_local);
      mpcx_si_sub (tmp, 2, tmp);
      mpcx_smul (h_local, h_local, tmp, exp_local);
   }

   if (overlap)
      mpcx_clear (h);
   mpcx_clear (tmp);
   mpcx_mv (h, h_local);
}