/*
test file for mul

Copyright (C) 2009, 2010, 2013 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 <stdio.h>
#include "mpfrcx.h"


static void error (mpcx_t h1, mpcx_t h2, mpcx_t f, mpcx_t g) {
   fprintf (stderr, "Error in mul: f * g yields h1 and h2 with\nf: ");
   mpcx_out_str (stderr, 10, 0, f);
   fprintf (stderr, "\ng: ");
   mpcx_out_str (stderr, 10, 0, g);
   fprintf (stderr, "\nh1: ");
   mpcx_out_str (stderr, 10, 0, h1);
   fprintf (stderr, "\nh2: ");
   mpcx_out_str (stderr, 10, 0, h2);
   fprintf (stderr, "\n");
   exit (1);
}


static void check_mul (mpcx_t f, mpcx_t g) {
   /* checks whether multiplication in place works; f and g must have the
      same precision */
   mpcx_t f2, g2, h;

   mpcx_init_set (f2, f);
   mpcx_init_set (g2, g);
   mpcx_init (h, 10, mpcx_get_prec (f));

   mpcx_mul (h, f, g);

   /* check for difference between normal and in place operation */
   mpcx_mul (f2, f2, g);
   if (mpcx_cmp (h, f2))
      error (h, f2, f, g);

   /* check for difference between normal and in place operation */
   mpcx_mul (g2, f, g2);
   if (mpcx_cmp (h, g2))
      error (h, g2, f, g);

   /* check whether f*g = g*f; should hold in floating point at least
      when using the fft */
   mpcx_mul (f2, g, f);
   if (mpcx_cmp (h, f2))
      error (h, f2, f, g);

   mpcx_clear (f2);
   mpcx_clear (g2);
   mpcx_clear (h);
}


static void check_mul_random (gmp_randstate_t state) {
   int deg;
   mpcx_t f, g;

   mpcx_init (f, 10, 103);
   mpcx_init (g, 11, 103);

   for (deg = -1; deg <= 1000; deg += 200) {
      mpcx_urandom (f, deg, state);
      mpcx_urandom (g, deg, state);
      check_mul (f, g);
      check_mul (f, f);
   }

   mpcx_clear (f);
   mpcx_clear (g);
}


static void check_unbalanced_mul () {
   const int degf = 4000, degg = 1000;
   const mpfr_prec_t prec = 200;
   mpcx_t f, g, h;
   mpfr_t coeff;
   unsigned long int coeff_int;
   int i;

   mpcx_init (f, degf + 1, prec);
   mpcx_init (g, degg + 1, prec);
   mpcx_init (h, degf + degg + 1, prec);
   mpfr_init2 (coeff, prec);

   mpcx_set_deg (f, degf);
   mpcx_set_deg (g, degg);
   for (i = 0; i <= degf; i++)
      mpc_set_ui (mpcx_get_coeff (f, i), i, MPC_RNDNN);
   for (i = 0; i <= degg; i++)
      mpc_set_ui (mpcx_get_coeff (g, i), degg + 1 - i, MPC_RNDNN);
   mpcx_mul (h, f, g);
   /* largest coefficient is that of x^4000, 1838837000 */
   mpc_abs (coeff, mpcx_get_coeff (h, degf), MPC_RNDNN);
   coeff_int = mpfr_get_ui (coeff, GMP_RNDN);

   if (coeff_int != 1838837000)
      exit (1);

   mpcx_clear (f);
   mpcx_clear (g);
   mpcx_clear (h);
   mpfr_clear (coeff);
}


int main (void) {
   gmp_randstate_t state;

   gmp_randinit_default (state);

   check_mul_random (state);

   gmp_randclear (state);

   check_unbalanced_mul ();

   return 0;
}