/*
mpcx_tower

functions working with towers of abelian extensions of \Q (\sqrt D),
where D is a negative discriminant

Copyright (C) 2017, 2018 Andreas Enge
Copyright (C) 2017 Jared Asuncion (INRIA)

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 int mpcx_is_galois_generator (mpc_t *embeddings, int n);
static void mpcx_tower_decomposition_rec (int levels, int *d,
   mpcx_t **W, mpcx_ptr V, mpc_t *roots);

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

mpfr_prec_t mpfr_coinciding_bits (mpfr_srcptr op1, mpfr_srcptr op2)
   /* The idea of this function is to return the number of leading bits in
      which op1 and op2 coincide.
      More precisely, if op1 and op2 are equal, then their minimal precision
      is returned.
      Otherwise if they have different signs, 0 is returned.
      Otherwise their difference is related to the larger of the two input
      exponents; this means that on input 1.000 and 0.111 in base 2, a
      value of 3 is returned.
      Moreover, the return value is capped at the minimal precision; this
      means that on input 1.000 and 1.0000000000001 a value of 4 is
      returned. */

{
   mpfr_t diff;
   mpfr_prec_t res, min_prec;

   if (!mpfr_number_p (op1) || !mpfr_number_p (op2))
      return (0);

   min_prec = MIN (mpfr_get_prec (op1), mpfr_get_prec (op2));
   mpfr_init2 (diff, 2);
   mpfr_sub (diff, op1, op2, MPFR_RNDU);
   if (mpfr_zero_p (diff))
      res = min_prec;
   else if (mpfr_zero_p (op1) || mpfr_zero_p (op2))
      res = 0;
   else if (   ( mpfr_signbit (op1) && !mpfr_signbit (op2))
            || (!mpfr_signbit (op1) &&  mpfr_signbit (op2)))
      res = 0;
   else
      res = MIN (min_prec,
                 MAX (mpfr_get_exp (op1), mpfr_get_exp (op2))
                    - mpfr_get_exp (diff));

   mpfr_clear (diff);

   return (res);
}

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

mpfr_prec_t mpc_coinciding_bits (mpc_srcptr op1, mpc_srcptr op2)
   /* Return the mininum of the coinciding bits of the real and imaginary
      parts of op1 and op2. */

{
   return (MIN (mpfr_coinciding_bits (mpc_realref (op1),
                                      mpc_realref (op2)),
                mpfr_coinciding_bits (mpc_imagref (op1),
                                      mpc_imagref (op2))));
}

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

static int mpcx_is_galois_generator (mpc_t *embeddings, int n)
   /* Given the n complex embeddings of an algebraic number in a Galois
      number field, the function returns 1 if the number is (modulo
      rounding errors) a generator of the field and 0 otherwise. */
{
   int i;
   mpfr_prec_t prec, min_prec;

   /* Every complex number in embeddings occurs the same number of times;
      so it is enough to check if embeddings [0] occurs a second time. */
   prec = mpc_get_prec (embeddings [0]);
   min_prec = prec / 2;
   for (i = 1; i < n; i++)
      /* We consider two numbers as equal if they coincide in at least half
         of their bits; since in the algorithms we have a large choice of
         generators, it is better to err by declaring distinct numbers as
         equal than the other way around. */
      if (mpc_coinciding_bits (embeddings [0], embeddings [1]) > min_prec)
         return (0);

   return (1);
}

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

void mpcx_tower_init (mpcx_tower_ptr twr, int levels, int* d,
   mpfr_prec_t prec)
   /* Initialise the tower field structure twr with levels field
      extensions, the relative degrees of which are passed via d.
      The common precision is given by prec. */
{
   int i, j, deg;

   twr->levels = levels;
   twr->d = (int *) malloc (levels * sizeof (int));
   deg = 1;
   for (i = 0; i < levels; i++) {
      twr->d [i] = d [i];
      deg *= d [i];
   }
   twr->deg = deg;
   twr->W = (mpcx_t **) malloc (levels * sizeof (mpcx_t *));
   deg = 1;
   for (i = 1; i < levels; i++) {
      twr->W [i] = (mpcx_t *) malloc ((d [i] + 1) * sizeof (mpcx_t));
      deg *= d [i - 1]; /* absolute degree of coefficients of W [i] */
      for (j = 0; j <= d [i]; j++)
         mpcx_init (twr->W [i][j], deg, prec);
   }
   twr->W [0] = (mpcx_t *) malloc (1 * sizeof (mpcx_t));
   mpcx_init (twr->W [0][0], d [0] + 1, prec);
}

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

void mpcx_tower_clear (mpcx_tower_ptr twr)
   /* Free all dynamically allocated components of a field tower
      structure. */
{
   int i, j;

   for (i = 1; i < twr->levels; i++) {
      for (j = 0; j <= twr->d [i]; j++)
         mpcx_clear (twr->W [i][j]);
      free (twr->W [i]);
   }
   free (twr->W [0]);
   free (twr->W);
   free (twr->d);
}
   
/**************************************************************************/

static void mpcx_tower_decomposition_rec (int levels, int *d,
   mpcx_t **W, mpcx_ptr V, mpc_t *roots)
   /* Compute recursively the tower field with levels >= 2 levels and
      degree sequence given by d, where the conjugates of a generating
      element x are given in the following order in roots: Let M / K be
      the total Galois extension of degree deg=m*n (the length of the
      array roots), and let L be the fixed field of the subgroup H of order
      m = d [levels-1], with [M:L]=m and [L:K]=n. Then it is assumed that
      the first m elements of roots are the conjugates of x under H, and
      more generally that roots is partitioned into n/m consecutive cosets
      of size m of xH by G/H. Moreover, these cosets are ordered
      recursively to be compatible with the normal series of the Galois
      group: The first few cosets correspond to the subgroup of G/H
      that determines the second relative extension from the top, and so on.
      The function determines a generator y of extension levels-2 (as an
      elementary symmetric function of the conjugates of x under H, for
      instance, if possible, its relative trace). It computes its absolute
      minimal polynomial V and fills W [levels-1] with the Hecke
      representation of the minimal polynomial of x over this extension.
      Then it fills roots with the n conjugates of y (respecting the order
      imposed by the normal series) and goes down the recursion.
      In the end, W [1], ..., W [levels-1] are filled with polynomials
      defining the levels-1 relative extensions from the top and V contains
      a polynomial defining the first absolute extension.
      All precisions should be the same; practically we take it from V.
   */
{
   int i, j, m, n, found;
   mpfr_prec_t prec;
   mpcx_t *U;
   mpc_t **vals;

   /* Set the precision and the degrees. */
   prec = mpcx_get_prec (V);
   m = d [levels-1];
   n = 1;
   for (i = 0; i < levels-1; i++)
      n *= d [i];

   /* Compute as floating point polynomials the n relative minimal
      polynomials U [k] of the conjugates of x over L. */
   U = (mpcx_t *) malloc (n * sizeof (mpcx_t));
   for (j = 0; j < n; j++) {
      mpcx_init (U [j], m+1, prec);
      mpcx_reconstruct_from_roots (U [j], roots + m*j, m);
   }

   /* y is one of the coefficients of U [0], and its conjugates are the
      corresponding coefficients of all the U [j]. Find a generator of
      the intermediate field among them, starting with the trace, which
      is the smallest candidate.
      To prepare the recursion, we replace the content of the array root
      by the conjugates of y. */
   found = 0;
   j = m-1;
   while (!found && j >= 0) {
      for (i = 0; i < n; i++)
         mpc_set (roots [i], mpcx_get_coeff (U [i], j), MPC_RNDNN);
      if (mpcx_is_galois_generator (roots, n))
         found = 1;
      else
         j--;
   }
   if (!found) {
      printf ("*** Houston, we have a problem!\n");
      printf ("*** No generator found in mpcx_tower_compute_rec.\n");
      printf ("*** Go back programming!\n");
      exit (1);
   }

   /* Compute V, the minimal polynomial of y, in a subproduct tree and
      simultaneously the W [levels-1], containing the Hecke representations
      of the coefficients of U [0]. Use W [levels-1] as a temporary
      variable and adjust later, since there is a shift by 1 in the result
      of mpcx_product_and_hecke_from_roots compared to the layout of our
      variables. */
   vals = (mpc_t **) malloc ((m+1) * sizeof (mpc_t *));
   for (j = 0; j <= m; j++) {
      vals [j] = (mpc_t *) malloc (n * sizeof (mpc_t));
      for (i = 0; i < n; i++)
         mpc_init2 (vals [j][i], prec);
   }
   for (i = 0; i < n; i++)
      mpc_set (vals [0][i], roots [i], MPC_RNDNN);
   for (j = 0; j < m; j++)
      for (i = 0; i < n; i++)
         mpc_set (vals [j+1][i], mpcx_get_coeff (U [i], j), MPC_RNDNN);
   for (j = 0; j < n; j++)
      mpcx_clear (U [j]);
   free (U);

   mpcx_product_and_hecke_from_roots (W [levels-1], vals, m+1, n);

   /* Shift the results back by 1 and replace W [levels-1][m] by the
      derivative of V. */
   mpcx_set (V, W [levels-1][0]);
   for (j = 0; j < m; j++)
      mpcx_set (W [levels-1][j], W [levels-1][j+1]);
   mpcx_derive (W [levels-1][m], V);

   for (j = 0; j <= m; j++) {
      for (i = 0; i < n; i++)
         mpc_clear (vals [j][i]);
      free (vals [j]);
   }
   free (vals);

   /* Recurse. */
   if (levels >= 3)
      mpcx_tower_decomposition_rec (levels - 1, d, W, V, roots);
}

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

void mpcx_tower_decomposition (mpcx_tower_ptr twr, mpc_t *roots)
   /* twr is an initialised structure for a Galois field tower, and roots
      a list of conjugates of a generating element in an order compatible
      with the normal series of the Galois group as explained in
      mpcx_tower_decomposition_rec; their precisions should be the same as
      that of twr.
      Compute the absolute and relative field representations in twr. */
{
   if (twr->levels == 1)
      mpcx_reconstruct_from_roots (twr->W [0][0], roots, twr->d [0]);
   else
      mpcx_tower_decomposition_rec (twr->levels, twr->d,
         twr->W, twr->W [0][0], roots);
}

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