/*
mpcx_treefromroots

high-level functions working with trees of polynomials, the leaves of
which are linear polynomials obtained from their roots

Copyright (C) 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-impl.h"

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

void mpcx_init_set_linear (mpcx_ptr h, mpc_srcptr z)
   /* Initialise h to the precision of z and set it to X-z. */

{
   mpfr_prec_t prec = mpc_get_prec (z);

   mpcx_init (h, 2, prec);
   mpcx_set_deg (h, 1);
   mpc_set_ui (h->coeff [1], 1, MPC_RNDNN);
   mpc_neg (h->coeff [0], z, MPC_RNDNN);
}

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

void mpcx_reconstruct_from_roots (mpcx_ptr rop, mpc_t *roots, int no)
   /* Takes an array of no elements and computes the polynomial having
      these as roots, that is, prod_{i=1}^n (X - roots[i]); the real
      and imaginary part of each element of roots must have the same
      precision. */

{
   int i;
   mpcx_t *fac;

   fac = (mpcx_t *) malloc (no * sizeof (mpcx_t));
   for (i = 0; i < no; i++)
      mpcx_init_set_linear (fac [i], roots [i]);

   mpcx_reconstruct (rop, fac, no);

   for (i = 0; i < no; i++)
      mpcx_clear (fac [i]);
   free (fac);
}

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

void mpcx_subproducttree_from_roots (mpcx_tree_ptr rop, mpc_t *roots,
   int no)
   /* Takes an array of no elements and initialises and computes the
      subproduct tree having the linear polynomials X - roots [i]
      as leaves. All roots are assumed to be of the same precision,
      which is used for the initialisation. */

{
   int i;
   mpcx_t *fac;
   mpfr_prec_t prec = mpc_get_prec (roots [0]);

   mpcx_tree_init (rop, no, prec);
   fac = (mpcx_t *) malloc (no * sizeof (mpcx_t));
   for (i = 0; i < no; i++)
      mpcx_init_set_linear (fac [i], roots [i]);

   mpcx_subproducttree (rop, fac);

   for (i = 0; i < no; i++)
      mpcx_clear (fac [i]);
   free (fac);
}

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

void mpcx_hecke_from_roots (mpcx_ptr rop, mpcx_tree_srcptr subproducts,
   mpc_t *vals)
   /* The function works exactly as mpcx_hecke, but assumes that the
      subproduct tree has been generated (for instance by a call to
      mpcx_subproducttree_from_roots) from monic linear polynomials;
      then vals is an array of constant polynomials, which may simply
      be passed as complex numbers. */

{
   int i, n = subproducts->no;
   mpfr_prec_t prec = mpc_get_prec (vals [0]);
   mpcx_t *fac;

   fac = (mpcx_t *) malloc (n * sizeof (mpcx_t));
   for (i = 0; i < n; i++) {
      mpcx_init (fac [i], 1, prec);
      mpcx_set_deg (fac [i], 0);
      mpcx_set_coeff (fac [i], 0, vals [i]);
   }

   mpcx_hecke (rop, subproducts, fac);

   for (i = 0; i < n; i++)
      mpcx_clear (fac [i]);
   free (fac);
}

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

void mpcx_product_and_hecke_from_roots (mpcx_t *rop, mpc_t **vals,
   int no_pols, int no_factors)
   /* The entries of val are expected to have the same precision. */

{
   int i, j;
   mpfr_prec_t prec = mpc_get_prec (vals [0][0]);
   mpcx_t **valx;

   valx = (mpcx_t **) malloc (no_pols * sizeof (mpcx_t *));
   for (i = 0; i < no_pols; i++)
      valx [i] = (mpcx_t *) malloc (no_factors * sizeof (mpcx_t));
   for (j = 0; j < no_factors; j++)
      mpcx_init_set_linear (valx [0][j], vals [0][j]);
   for (i = 1; i < no_pols; i++)
      for (j = 0; j < no_factors; j++) {
         mpcx_init (valx [i][j], 1, prec);
         mpcx_set_deg (valx [i][j], 0);
         mpcx_set_coeff (valx [i][j], 0, vals [i][j]);
      }

   mpcx_product_and_hecke (rop, valx, no_pols, no_factors);

   for (i = 0; i < no_pols; i++) {
      for (j = 0; j < no_factors; j++)
         mpcx_clear (valx [i][j]);
      free (valx [i]);
   }
   free (valx);
}

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