flint-sys 0.9.0

/*
    Copyright (C) 2011 Andy Novocin
    Copyright (C) 2011, 2016 William Hart

    This file is part of FLINT.

    FLINT is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
*/

#include <math.h>
#include "longlong.h"
#include "fmpz.h"
#include "fmpz_vec.h"
#include "fmpz_mat.h"
#include "fmpz_poly.h"
#include "fmpz_poly_factor.h"
#include "fmpz_lll.h"

static slong _heuristic_van_hoeij_starting_precision(const fmpz_poly_t f,
                                                            slong r, ulong p)
{
   fmpz_t lead_b, trail_b;
   slong min_b, a;

   fmpz_init(lead_b);
   fmpz_init(trail_b);

   fmpz_poly_CLD_bound(lead_b, f, f->length - 2);
   fmpz_poly_CLD_bound(trail_b, f, 0);

   min_b = FLINT_MIN(fmpz_bits(lead_b), fmpz_bits(trail_b));

   a = (slong) ((2.5*r + min_b)*log(2) + log(f->length)/2.0)/log(p);

   fmpz_clear(trail_b);
   fmpz_clear(lead_b);

   return a;
}

/*
   resize the matrix M to num_rows, where num_rows <= M->r
*/
/* todo: proper inplace realloc */
static void fmpz_mat_van_hoeij_resize_matrix(fmpz_mat_t M, slong num_rows)
{
    slong i;
    fmpz_mat_t T;

    FLINT_ASSERT(num_rows <= M->r);

    fmpz_mat_init(T, num_rows, M->c);
    for (i = 0; i < num_rows; i++)
        _fmpz_vec_swap(fmpz_mat_row(T, i), fmpz_mat_row(M, i), M->c);
    fmpz_mat_swap(M, T);
    fmpz_mat_clear(T);
}

void fmpz_poly_factor_van_hoeij(fmpz_poly_factor_t final_fac,
         const nmod_poly_factor_t fac, const fmpz_poly_t f, slong exp, ulong p)
{
   fmpz_poly_factor_t lifted_fac;
   fmpz_mat_t M;
   fmpz_t fp, P, B, lc, bound_sum;
   slong i, r = fac->num;
   slong bit_r = FLINT_MAX(r, 20);
   slong U_exp, a, next_col, num_data_cols, num_coeffs;
   slong prev_num_coeffs, prev_exp, N, worst_exp, num_rows;
   ulong sqN;
   fmpz_poly_t * v, * w;
   fmpz_mat_t col, data;
   slong * link;
   int hensel_loops, do_lll;
   fmpz_lll_t fl;

   /* set to identity */
   fmpz_mat_init(M, r, r);
   for (i = 0; i < r; i++)
      fmpz_one(fmpz_mat_entry(M, i, i));

   /* we prescale the identity matrix by 2^U_exp */
   U_exp = (slong) FLINT_BIT_COUNT(bit_r);

   fmpz_mat_scalar_mul_2exp(M, M, U_exp);

   /* compute Mignotte bound */
   fmpz_init(B);

   fmpz_poly_factor_mignotte(B, f);
   /*
      bound adjustment, we multiply true factors (which might be
      monic) by the leading coefficient of f in the implementation
      below
   */
   fmpz_mul(B, B, f->coeffs + f->length - 1);
   fmpz_abs(B, B);
   fmpz_mul_ui(B, B, 2);
   fmpz_add_ui(B, B, 1);
   a = fmpz_clog_ui(B, p);

   /* compute heuristic starting precision */
   a = FLINT_MIN(a, _heuristic_van_hoeij_starting_precision(f, r, p));

   /* start Hensel lift */
   fmpz_poly_factor_init(lifted_fac);

   v = flint_malloc((2*r - 2)*sizeof(fmpz_poly_t));
   w = flint_malloc((2*r - 2)*sizeof(fmpz_poly_t));
   link = flint_malloc((2*r - 2)*sizeof(slong));

   for (i = 0; i < 2*r - 2; i++)
   {
      fmpz_poly_init(v[i]);
      fmpz_poly_init(w[i]);
   }

   prev_exp = _fmpz_poly_hensel_start_lift(lifted_fac, link, v, w, f, fac, a);

   /* compute bound */
   fmpz_set_ui(B, r + 1);
   fmpz_mul_2exp(B, B, 2*U_exp);

   /* compute leading coefficient */
   N = f->length - 1;
   sqN = (ulong) sqrt(N);
   fmpz_init(lc);
   fmpz_set(lc, f->coeffs + N);

   /* main hensel loop */
   hensel_loops = 0;
   fmpz_init(P);
   fmpz_init(fp);
   fmpz_set_ui(fp, p);
   fmpz_pow_ui(P, fp, a);

   fmpz_init(bound_sum);
   fmpz_mat_init(col, r, 1);
   fmpz_lll_context_init_default(fl);

   while (!fmpz_poly_factor_van_hoeij_check_if_solved(M, final_fac, lifted_fac, f, P, exp, lc))
   {
      if (hensel_loops < 3 && 3*r > N + 1)
         num_coeffs = r > 200 ? 50 : 30;
      else
         num_coeffs = 10;

      num_coeffs = FLINT_MIN(num_coeffs, (N + 1)/2);
      prev_num_coeffs = 0;

      do {
         fmpz_mat_init(data, r + 1, 2*num_coeffs);
         num_data_cols = _fmpz_poly_factor_CLD_mat(data, f, lifted_fac, P, num_coeffs);

         for (next_col = prev_num_coeffs; next_col < num_data_cols - prev_num_coeffs; next_col++)
         {
            /* we alternate taking columns from the right and left */
            slong alt_col, diff = next_col - prev_num_coeffs;

            if ((diff % 2) == 0)
               alt_col = prev_num_coeffs + diff/2;
            else
               alt_col = num_data_cols - prev_num_coeffs - (diff + 1)/2;

            fmpz_mul_ui(bound_sum, fmpz_mat_entry(data, r, alt_col), sqN);
            worst_exp = fmpz_bits(bound_sum);

            for (i = 0; i < r; i++)
               fmpz_set(fmpz_mat_entry(col, i, 0), fmpz_mat_entry(data, i, alt_col));

            do_lll = fmpz_mat_next_col_van_hoeij(M, P, col, worst_exp, U_exp);

            if (do_lll)
            {
                /* flint_printf("LLL %wd %wd  %wd\n", M->r, M->c, fmpz_mat_max_bits(M)); */
               num_rows = fmpz_lll_wrapper_with_removal_knapsack(M, NULL, B, fl);

               fmpz_mat_van_hoeij_resize_matrix(M, num_rows);

               if (fmpz_poly_factor_van_hoeij_check_if_solved(M, final_fac, lifted_fac, f, P, exp, lc))
               {
                  fmpz_mat_clear(data);
                  goto cleanup;
               }
            }
         }

         prev_num_coeffs = num_coeffs;
         num_coeffs = FLINT_MIN(2*num_coeffs, (N + 1)/2);
         fmpz_mat_clear(data);

      } while (num_coeffs != prev_num_coeffs);

      hensel_loops++;

      prev_exp = _fmpz_poly_hensel_continue_lift(lifted_fac, link, v, w, f, prev_exp, a, 2*a, fp);

      a = 2*a;
      fmpz_pow_ui(P, fp, a);
   }

cleanup:

   fmpz_clear(lc);
   fmpz_clear(fp);
   fmpz_clear(P);
   fmpz_clear(B);
   fmpz_mat_clear(col);
   fmpz_mat_clear(M);
   fmpz_clear(bound_sum);
   fmpz_poly_factor_clear(lifted_fac);

   for (i = 0; i < 2*r - 2; i++)
   {
      fmpz_poly_clear(v[i]);
      fmpz_poly_clear(w[i]);
   }
   flint_free(v);
   flint_free(w);
   flint_free(link);
}