/*
mpfrx_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 mpfrx_rem_naive (mpfrx_ptr r, mpfrx_srcptr f, mpfrx_srcptr g);
static void mpfrx_rem_newton (mpfrx_ptr r, mpfrx_srcptr f, mpfrx_srcptr g);
static void mpfrx_rev (mpfrx_ptr h, mpfrx_srcptr f, const int exp);
static void mpfrx_smul (mpfrx_ptr h, mpfrx_ptr f, mpfrx_ptr g, const int exp);
static void mpfrx_sinv (mpfrx_ptr h, mpfrx_ptr f, const int exp);

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

void mpfrx_rem (mpfrx_ptr r, mpfrx_srcptr f, mpfrx_srcptr g)

{
   if (f->deg < g->deg)
      mpfrx_set (r, f);
   else
      if (g->deg >= 32 && f->deg >= 63)
         mpfrx_rem_newton (r, f, g);
      else
         mpfrx_rem_naive (r, f, g);
}

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

static void mpfrx_rem_naive (mpfrx_ptr r, mpfrx_srcptr f, mpfrx_srcptr g)
   /* computes the remainder of f divided by g by the naive algorithm for    */
   /* deg f >= deg g                                                         */
{
   mpfr_t tmp;
   int i, j;

   if (mpfr_cmp_si (g->coeff [g->deg], 1) == 0) {
      mpfr_init2 (tmp, mpfr_get_prec (f->coeff [0]));
      /* the only case we actually use for multievaluation */
      mpfrx_set (r, f);
      for (i = f->deg - g->deg; i >= 0; i--)
         for (j = 0; j < g->deg; j++) {
            mpfr_mul (tmp, r->coeff [g->deg + i], g->coeff [j], GMP_RNDN);
            mpfr_sub (r->coeff [i + j], r->coeff [i + j], tmp, GMP_RNDN);
         }
      r->deg = g->deg - 1;
      /* To save memory, one might want to do a realloc on r here:
      mpfrx_realloc (r, r->deg + 1);
      */
      mpfr_clear (tmp);
   }
   else {
      printf ("*** Houston, we have a problem!\n");
      printf ("*** Calling mpfrx_rem_naive with a non-monic divisor.\n");
      printf ("*** Go back programming!\n");
      exit (1);
   }
}

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

static void mpfrx_rem_newton (mpfrx_ptr r, mpfrx_srcptr f, mpfrx_srcptr g) {
   /* computes the remainder of f divided by g by Newton iterations for      */
   /* deg f >= deg g                                                         */

   mpfrx_t q, tmp;

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

   mpfrx_rev (tmp, f, 0);
   mpfrx_rev (q, g, 0);
   mpfrx_sinv (q, q, f->deg - g->deg + 1);
   mpfrx_smul (q, tmp, q, f->deg - g->deg + 1);
   mpfrx_rev (q, q, f->deg - g->deg);

   mpfrx_mul (tmp, q, g);
   mpfrx_sub (r, f, tmp);
   if (r->deg >= g->deg)
      r->deg = g->deg - 1;

   mpfrx_clear (q);
   mpfrx_clear (tmp);
}

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

static void mpfrx_rev (mpfrx_ptr h, mpfrx_srcptr f, const int exp) {
   int    overlap, exp_local, i;
   mpfrx_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 'mpfrx_rev'.\n");
      exit (1);
   }

   overlap = (h == f);
   if (overlap)
      mpfrx_init (h_local, exp_local + 1, h->prec);
   else {
      mpfrx_mv (h_local, h);
      if (h_local->size < exp_local + 1)
         mpfrx_realloc (h_local, exp_local + 1);
   }

   h_local->deg = exp_local;
   for (i = 0; i <= f->deg; i++)
      mpfr_set (h_local->coeff [exp_local - i], f->coeff [i], GMP_RNDN);
   for (i = exp_local - f->deg - 1; i >= 0; i--)
      mpfr_set_ui (h_local->coeff [i], 0, GMP_RNDN);

   while (h_local->deg >= 0
         && mpfr_cmp_si (h_local->coeff [h_local->deg], 0) == 0)
      h_local->deg--;
   if (overlap)
      mpfrx_clear (h);
   mpfrx_mv (h, h_local);
}

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

static void mpfrx_smul (mpfrx_ptr h, mpfrx_ptr f, mpfrx_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;
   mpfrx_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 && mpfr_cmp_si (h->coeff [h->deg], 0) == 0)
         h->deg--;
   }
}

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

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

   int    overlap;
   mpfrx_t h_local, tmp;
   int    exp_local;

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

   exp_local = 1;
   h_local->deg = 0;
   mpfr_ui_div (h_local->coeff [0], 1, f->coeff [0], GMP_RNDN);
   while (exp_local < exp) {
      exp_local = (2 * exp_local < exp ? 2 * exp_local : exp);
      mpfrx_smul (tmp, h_local, f, exp_local);
      mpfrx_si_sub (tmp, 2, tmp);
      mpfrx_smul (h_local, h_local, tmp, exp_local);
   }

   if (overlap)
      mpfrx_clear (h);
   mpfrx_clear (tmp);
   mpfrx_mv (h, h_local);
}