/*
test file for tree functions

Copyright (C) 2017, 2018 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.h"

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

static int real_get_si (long int *z, mpfr_srcptr f)
   /* Round the real to an integer value; the return value reflects the
      success of the operation, in the sense that the difference is
      sufficiently small. */

{
   mpfr_t diff;
   mpfr_prec_t prec;
   mpfr_exp_t e;
   int ok;

   prec = mpfr_get_prec (f);
   e = mpfr_get_exp (f);
   if (prec <= e)
      /* The number is trivially an integer, but with no fractional part,
         so the rounding error has become invisible. */
      return (0);

   mpfr_init2 (diff, 2);
   *z = mpfr_get_si (f, MPFR_RNDN);
   mpfr_sub_si (diff, f, *z, MPFR_RNDN);

   ok = (mpfr_sgn (diff) == 0 || mpfr_get_exp (diff) <= -10);
      /* We arbitrarily accept a rounding error of at most 2^(-10). */

   mpfr_clear (diff);

   return (ok);
}

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

static int example1 (void) {

   int i, j, k;
   mpfr_prec_t prec = 128;
   mpc_t *roots;
   int deg = 20, d [3] = { 5, 2, 2 };
   int ok;
   mpcx_tower_t twr;
   mpc_ptr coeff;
   long int z;
   const long int W [3][3][10]
      = {{{ -73, -80, -74, -39, -6, 1, 0, 0, 0, 0},
          { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
          { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
         {{ -140, 53, 54, 12, 1, 0, 0, 0, 0, 0},
          { 365, 320, 222, 78, 6, 0, 0, 0, 0, 0},
          { -80, -148, -117, -24, 5, 0, 0, 0, 0, 0}},
         {{ -95, 55, 228, -250, -47, 160, 20, -142, 80, -12},
          { -110, -279, 1224, -1414, 624, -60, 116, -186, 76, -6},
          { 31, -306, 606, -416, 60, -174, 434, -304, 54, 10 }}};

   roots = (mpc_t *) malloc (deg * sizeof (mpc_t));
   for (i = 0; i < deg; i++)
      mpc_init2 (roots [i], prec);

   mpc_set_str (roots [ 0], "11.537612002601293729555643926712700128", 10, MPC_RNDNN);
   mpc_set_str (roots [ 1], "-1.2359811407720793940569330632262337764", 10, MPC_RNDNN);
   mpc_set_str (roots [ 2], "-0.88158412095480496703824738845507654109", 10, MPC_RNDNN);
   mpc_set_str (roots [ 3], "1.0537757586472475214611033030981796170", 10, MPC_RNDNN);
   mpc_set_str (roots [ 4], "(-1.3738815073324657798749662481447867991 -1.1528234664875328372707537955648321806)", 10, MPC_RNDNN);
   mpc_set_str (roots [ 5], "(0.15539825039129532339792538912105083438 0.83255682211942350580400245155255066129)", 10, MPC_RNDNN);
   mpc_set_str (roots [ 6], "(0.23751391794267422186896086399600757039 -0.66425231309624577378237568506293771302)", 10, MPC_RNDNN);
   mpc_set_str (roots [ 7], "(0.75323765599341107222341710074633201727 -0.24354073806061626168053242233604424063)", 10, MPC_RNDNN);
   mpc_set_str (roots [ 8], "(-0.66231668017784407796844395136033368855 0.70148728032360712487311562970674213590)", 10, MPC_RNDNN);
   mpc_set_str (roots [ 9], "(-0.62514835842316383925952894202832472183 -0.49160598090266240182413906785025006010)", 10, MPC_RNDNN);
   mpc_set_str (roots [10], "(-0.22655403352361409872463829988445090319 -0.81762769695565820067029579086297358680)", 10, MPC_RNDNN);
   mpc_set_str (roots [11], "(-0.49516049463112126662350930151027902325 -0.048769033268922442525849010844560810424)", 10, MPC_RNDNN);
   mpc_conj (roots [12], roots [10], MPC_RNDNN);
   mpc_conj (roots [13], roots [11], MPC_RNDNN);
   mpc_conj (roots [14], roots [ 8], MPC_RNDNN);
   mpc_conj (roots [15], roots [ 9], MPC_RNDNN);
   mpc_conj (roots [16], roots [ 6], MPC_RNDNN);
   mpc_conj (roots [17], roots [ 7], MPC_RNDNN);
   mpc_conj (roots [18], roots [ 4], MPC_RNDNN);
   mpc_conj (roots [19], roots [ 5], MPC_RNDNN);

   mpcx_tower_init (twr, 3, d, prec);
   mpcx_tower_decomposition (twr, roots);

   /* Check the degrees. */
   ok = mpfrx_get_deg (twr->W [0][0]) == 5;
   for (j = 0; j < 3; j++)
      ok = ok && mpfrx_get_deg (twr->W [1][j]) == 4;
   for (j = 0; j < 3; j++)
      ok = ok && mpfrx_get_deg (twr->W [2][j]) == 9;

   /* Check whether the result is an integral polynomial without
      imaginary part and with the correct real part. */
   for (i = 0; i < 3; i++)
      for (j = 0; j <= (i == 0 ? 0 : d [i]); j++)
         for (k = 0; k <= mpcx_get_deg (twr->W [i][j]); k++) {
            coeff = mpcx_get_coeff (twr->W [i][j], k);
            ok = ok && real_get_si (&z, mpc_realref (coeff))
                    && z == W [i][j][k]
                    && real_get_si (&z, mpc_imagref (coeff))
                    && z == 0;
      }
   
   mpcx_tower_clear (twr);

   for (i = 0; i < deg; i++)
      mpc_clear (roots [i]);
   free (roots);

   return (ok);
}

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

static int example2 (void) {
   /* only one level */

   int i, k;
   mpfr_prec_t prec = 128;
   mpc_t *roots;
   int deg = 3, d [1] = { 3 };
   int ok;
   mpcx_tower_t twr;
   mpc_ptr coeff;
   long int z;
   const long int W [4] = { -2, 0, 0, 1 };

   roots = (mpc_t *) malloc (deg * sizeof (mpc_t));
   for (i = 0; i < deg; i++)
      mpc_init2 (roots [i], prec);

   mpc_set_str (roots [ 0], "1.2599210498948731647672106072782283506", 10, MPC_RNDNN);
   mpc_set_str (roots [ 1], "(-0.62996052494743658238360530363911417528 1.0911236359717214035600726141898088813)", 10, MPC_RNDNN);
   mpc_conj (roots [2], roots [1], MPC_RNDNN);

   mpcx_tower_init (twr, 1, d, prec);
   mpcx_tower_decomposition (twr, roots);

   /* Check whether the result is an integral polynomial without
      imaginary part and with the correct real part. */
   ok = mpcx_get_deg (twr->W [0][0]) == 3;
   for (k = 0; k < 4; k++) {
      coeff = mpcx_get_coeff (twr->W [0][0], k);
      ok = ok && real_get_si (&z, mpc_realref (coeff))
         && z == W [k]
         && real_get_si (&z, mpc_imagref (coeff))
         && z == 0;
   }
   
   mpcx_tower_clear (twr);

   for (i = 0; i < deg; i++)
      mpc_clear (roots [i]);
   free (roots);

   return (ok);
}

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

int main (void) {
   return !(example1 () && example2 ());
}