/*
test file for tree functions

Copyright (C) 2011, 2012 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 check_reconstruct_multieval () {

   const int n=12;
   const mpfr_prec_t prec = 100;
   int i;
   mpc_t *values, *args;
   mpcx_t f;

   values = (mpc_t *) malloc (n * sizeof (mpc_t));
   args = (mpc_t *) malloc (n * sizeof (mpc_t));
   for (i = 0; i < n; i++) {
      mpc_init2 (values [i], prec);
      mpc_init2 (args [i], prec);
      mpc_set_ui (args [i], (unsigned long int) i, MPC_RNDNN);
   }
   
   mpcx_init (f, n+1, prec);
   mpcx_reconstruct_from_roots (f, args, n);
   mpcx_multieval (values, args, n, f);
   mpcx_clear (f);

   for (i = 0; i < n; i++) {
      if (mpc_cmp (values [i], args [0]) != 0) {
         fprintf (stderr, "Error in reconstruct or multieval: Value should be 0, but is\n");
         mpc_out_str (stdout, 10, 0, values [i], MPC_RNDNN);
         fprintf (stderr, "\n");
         exit (1);
      }
      mpc_clear (values [i]);
      mpc_clear (args [i]);
   }
   free (values);
   free (args);
}


static void
check_hecke (void) {
   const mpfr_prec_t prec = 100;
   int i, j;
   mpcx_t factor [3], f;
   mpcx_tree_t subprods;
   mpcx_t val [3], res;
   mpcx_t **vals, *ress;

   mpcx_init (f, 5, prec);
   mpcx_init (res, 4, prec);
   for (j = 0; j < 3; j++) {
      mpcx_init (factor [j], 3, prec);
      mpcx_init (val [j], 2, prec);
   }
   ress = (mpcx_t *) malloc (3 * sizeof (mpcx_t));
   vals = (mpcx_t **) malloc (3 * sizeof (mpcx_t *));
   for (i = 0; i < 3; i++) {
      mpcx_init (ress [i], 5, prec);
      vals [i] = (mpcx_t *) malloc (3 * sizeof (mpcx_t));
      for (j = 0; j < 3; j++)
         mpcx_init (vals [i][j], 3, prec);
   }
   mpcx_tree_init (subprods, 3, prec);

   mpcx_set_deg (factor [0], 2);
   mpc_set_ui (factor [0]->coeff [2], 1ul, MPC_RNDNN);
   mpc_set_ui (factor [0]->coeff [1], 0ul, MPC_RNDNN);
   mpc_set_ui (factor [0]->coeff [0], 1ul, MPC_RNDNN);
   mpcx_set_deg (factor [1], 1);
   mpc_set_ui (factor [1]->coeff [1], 1ul, MPC_RNDNN);
   mpc_set_si (factor [1]->coeff [0], -1l, MPC_RNDNN);
   mpcx_set_deg (factor [2], 1);
   mpc_set_ui (factor [2]->coeff [1], 1ul, MPC_RNDNN);
   mpc_set_si (factor [2]->coeff [0], -2l, MPC_RNDNN);
   mpcx_set_deg (val [0], 1);
   mpc_set_ui (val [0]->coeff [1], 2ul, MPC_RNDNN);
   mpc_set_si (val [0]->coeff [0], -2l, MPC_RNDNN);
   mpcx_set_deg (val [1], 0);
   mpc_set_ui (val [1]->coeff [0], 3ul, MPC_RNDNN);
   mpcx_set_deg (val [2], 0);
   mpc_set_ui (val [2]->coeff [0], 4ul, MPC_RNDNN);
   mpcx_set_deg (res, 3);
   mpc_set_ui (res->coeff [3], 9ul, MPC_RNDNN);
   mpc_set_si (res->coeff [2], -18l, MPC_RNDNN);
   mpc_set_ui (res->coeff [1], 17ul, MPC_RNDNN);
   mpc_set_si (res->coeff [0], -14l, MPC_RNDNN);

   mpcx_subproducttree (subprods, factor);
   mpcx_hecke (f, subprods, val);

   if (mpcx_cmp (f, res) != 0) {
      fprintf (stderr, "Error in hecke\n");
      printf ("f   "); mpcx_out_str (stdout, 10, 0, f); printf ("\n");
      printf ("res "); mpcx_out_str (stdout, 10, 0, res); printf ("\n");
      exit (1);
   }

   for (j = 0; j < 3; j++)
      mpcx_set (vals [0][j], factor [j]);
   for (i = 1; i < 3; i++)
      for (j = 0; j < 3; j++)
         mpcx_set (vals [i][j], val [j]);

   mpcx_product_and_hecke (ress, vals, 3, 3);
   mpcx_tree_get_root (f, subprods);
   if (mpcx_cmp (f, ress [0]) != 0) {
      fprintf (stderr, "Error in product_and_hecke\n");
      printf ("f        "); mpcx_out_str (stdout, 10, 0, f); printf ("\n");
      printf ("ress [0] "); mpcx_out_str (stdout, 10, 0, ress [0]); printf ("\n");
      exit (1);
   }
   if (mpcx_cmp (res, ress [1]) != 0) {
      fprintf (stderr, "Error in product_and_hecke\n");
      printf ("res      "); mpcx_out_str (stdout, 10, 0, res); printf ("\n");
      printf ("ress [1] "); mpcx_out_str (stdout, 10, 0, ress [1]); printf ("\n");
      exit (1);
   }
   if (mpcx_cmp (res, ress [2]) != 0) {
      fprintf (stderr, "Error in product_and_hecke\n");
      printf ("res      "); mpcx_out_str (stdout, 10, 0, res); printf ("\n");
      printf ("ress [2] "); mpcx_out_str (stdout, 10, 0, ress [2]); printf ("\n");
      exit (1);
   }

   mpcx_clear (f);
   mpcx_clear (res);
   for (j = 0; j < 3; j++) {
      mpcx_clear (factor [j]);
      mpcx_clear (val [j]);
   }
   for (i = 0; i < 3; i++) {
      mpcx_clear (ress [i]);
      for (j = 0; j < 3; j++)
         mpcx_clear (vals [i][j]);
      free (vals [i]);
   }
   free (ress);
   free (vals);
   mpcx_tree_clear (subprods);
}


int main (void) {

   check_reconstruct_multieval ();
   check_hecke ();

   return 0;
}