flint-sys 0.9.0

/*
    Copyright (C) 2013 Sebastian Pancratz
    Copyright (C) 2015 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 "ulong_extras.h"
#include "nmod_mat.h"
#include "nmod_poly.h"
#include "fmpz.h"
#include "fmpz_vec.h"
#include "fmpz_poly.h"
#include "fmpz_mat.h"

static slong _fmpz_mat_minpoly_small(fmpz * rop, const fmpz_mat_t op)
{
    slong len = 0;

    if (op->r == 0)
    {
        fmpz_one(rop + 0);
        len = 1;
    }
    else if (op->r == 1)
    {
        fmpz_one(rop + 1);
        fmpz_neg(rop + 0, fmpz_mat_entry(op, 0, 0));
        len = 2;
    }

    return len;
}

static void _fmpz_mat_bound_ovals_of_cassini(fmpz_t b, const fmpz_mat_t op)
{
   slong n = op->r, i, j;
   fmpz * v1, * v2;
   fmpz_t t, q, r1, r2;

   fmpz_init(t);
   fmpz_init(q);
   fmpz_init(r1);
   fmpz_init(r2);

   v1 = _fmpz_vec_init(n);
   v2 = _fmpz_vec_init(n);

   /* |A| [1,1,...,1]^T */
   for (i = 0; i < n; i++)
   {
      for (j = 0; j < n; j++)
      {
         if (fmpz_sgn(fmpz_mat_entry(op, i, j)) >= 0)
            fmpz_add(v1 + i, v1 + i, fmpz_mat_entry(op, i, j));
         else
            fmpz_sub(v1 + i, v1 + i, fmpz_mat_entry(op, i, j));
      }
   }
   /* |A|^t * |A| * [1,1,...,1]^T */
   for (i = 0; i < n; i++)
   {
      for (j = 0; j < n; j++)
      {
         if (fmpz_sgn(fmpz_mat_entry(op, j, i)) >= 0)
            fmpz_addmul(v2 + i, v1 + j, fmpz_mat_entry(op, j, i));
         else
            fmpz_submul(v2 + i, v1 + j, fmpz_mat_entry(op, j, i));
      }
   }

   _fmpz_vec_clear(v1, n);
   v1 = v2;

   for (i = 0; i < n; i++)
   {
      fmpz_zero(t);
      /* q_i */
      for (j = 1; j < n; j++)
         fmpz_addmul(t, fmpz_mat_entry(op, i, j), fmpz_mat_entry(op, i, j));

      if (fmpz_cmp(t, q) > 0)
         fmpz_set(q, t);

      /* r_i */
      fmpz_sub(t, v1 + i, t);

      if (fmpz_cmp(t, r2) > 0)
      {
         fmpz_swap(t, r2);

         if (fmpz_cmp(r2, r1) > 0)
            fmpz_swap(r2, r1);
      }
   }

   fmpz_mul(r1, r1, r2);

   fmpz_sqrtrem(b, r2, r1);

   if (!fmpz_is_zero(r2))
      fmpz_add_ui(b, b, 1);

   fmpz_add(b, b, q);

   _fmpz_vec_clear(v1, n);
   fmpz_clear(r1);
   fmpz_clear(r2);
   fmpz_clear(t);
   fmpz_clear(q);
}

static inline double _log2e(const double x, slong exp)
{
    return log(x) * 1.44269504088896340736 + exp;
}

/*
    If $A$ is an $n \times n$ matrix with spectral radius
    bound by b, the coefficients of the minimal polynomial have
    at most $\ceil{d\log_2(b)}$ bits if $d \leq b$. Otherwise
    if $d > 0$ it has at most
    $\min{\ceil{d/2\log_2(bd)}, \ceil{d\log_2(2b)}}$
    bits, where $d$ is the degree of the minimal polynomial.
    See Lemma 3.1 in Dumas, "Bounds on the coefficients of the
    characteristic and minimal polynomials", 2007.
*/
static slong
_fmpz_mat_minpoly_bound_bits(const fmpz_mat_t op)
{
    const slong n = op->r;
    slong bound;

    fmpz_t b;
    fmpz_init(b);

    _fmpz_mat_bound_ovals_of_cassini(b, op);
    fmpz_add_ui(b, b, fmpz_is_zero(b));

    double bb, b1, b2;
    slong exp;

    bb = fmpz_get_d_2exp(&exp, b);

    fmpz_clear(b);

    if (fmpz_cmp_ui(b, n) >= 0)
    {
        b1 = _log2e(bb, exp);
    }
    else
    {
        b1 = _log2e(bb * n, exp) * 0.5;
        b2 = _log2e(bb * 2, exp);
        b1 = FLINT_MIN(b1, b2);
    }

    bound = ceil(n * b1 * (1 + 1e-10));

    return bound;
}

slong _fmpz_mat_minpoly_modular(fmpz * rop, const fmpz_mat_t op)
{
    const slong n = op->r;
    slong len = 0, oldlen = 0;

    if (n < 2)
    {
        return _fmpz_mat_minpoly_small(rop, op);
    }
    else
    {
        /*
            If $A$ is an $n \times n$ matrix with spectral radius
            bound by b, the coefficients of the minimal polynomial have
            at most $\ceil{d\log_2(b)}$ bits if $d \leq b$. Otherwise
            if $d > 0$ it has at most
            $\min{\ceil{d/2\log_2(bd)}, \ceil{d\log_2(2b)}}$
            bits, where $d$ is the degree of the minimal polynomial.
            See Lemma 3.1 in Dumas, "Bounds on the coefficients of the
            characteristic and minimal polynomials", 2007.
        */
        ulong bound;

        slong pbits  = FLINT_BITS - 1, i, j;
        ulong p = (UWORD(1) << pbits);
        ulong * P, * Q;

        fmpz_mat_t v1, v2, v3;
        fmpz * rold;
        fmpz_t m;

        if (fmpz_mat_is_zero(op))
        {
           fmpz_set_ui(rop + 0, 0);
           fmpz_set_ui(rop + 1, 1);
           return 2;
        }

        /* Determine the bound in bits */
        bound = _fmpz_mat_minpoly_bound_bits(op);
        /* Allow for signs */
        bound += 1;

        P = (ulong *) flint_calloc(n, sizeof(ulong));
        Q = (ulong *) flint_calloc(n, sizeof(ulong));
        rold = (fmpz *) _fmpz_vec_init(n + 1);
        fmpz_mat_init(v1, n, 1);
        fmpz_mat_init(v2, n, 1);
        fmpz_mat_init(v3, n, 1);

        fmpz_init_set_ui(m, 1);

        oldlen = 0;
        len = 0;

        for ( ; fmpz_bits(m) <= bound; )
        {
            nmod_mat_t mat;
            nmod_poly_t poly;

            p = n_nextprime(p, 0);

            nmod_mat_init(mat, n, n, p);
            nmod_poly_init(poly, p);

            for (i = 0; i < n; i++)
               P[i] = 0;

            fmpz_mat_get_nmod_mat(mat, op);
            nmod_mat_minpoly_with_gens(poly, mat, P);

            len = poly->length;

            if (oldlen != 0 && len > oldlen)
            {
               /* all previous primes were bad, discard */

               fmpz_one(m);
               oldlen = len;

               for (i = 0; i < n + 1; i++)
                  fmpz_zero(rop + i);

               for (i = 0; i < n; i++)
                  Q[i] = 0;
            } else if (len < oldlen)
            {
               /* this prime was bad, skip */

               nmod_mat_clear(mat);
               nmod_poly_clear(poly);

               continue;
            }

            for (i = 0; i < n; i++)
               Q[i] |= P[i];

            _fmpz_poly_CRT_ui(rop, rop, n + 1, m, poly->coeffs,
                              poly->length, poly->mod.n, poly->mod.ninv, 1);

            fmpz_mul_ui(m, m, p);

            /* check if stabilised */
            for (i = 0; i < len; i++)
            {
               if (!fmpz_equal(rop + i, rold + i))
                  break;
            }

            for (j = 0; j < len; j++)
               fmpz_set(rold + j, rop + j);

            if (i == len) /* stabilised */
            {
               for (i = 0; i < n; i++)
               {
                  if (Q[i] == 1)
                  {
                     fmpz_mat_zero(v1);
                     fmpz_mat_zero(v3);

                     fmpz_set_ui(fmpz_mat_entry(v1, i, 0), 1);

                     for (j = 0; j < len; j++)
                     {
                        fmpz_mat_scalar_mul_fmpz(v2, v1, rop + j);
                        fmpz_mat_add(v3, v3, v2);

                        if (j != len - 1)
                        {
                           fmpz_mat_mul(v2, op, v1);
                           fmpz_mat_swap(v1, v2);
                        }
                     }

                     /* check f(A)v = 0 */
                     for (j = 0; j < n; j++)
                     {
                        if (!fmpz_is_zero(fmpz_mat_entry(v3, j, 0)))
                            break;
                     }

                     if (j != n)
                        break;
                  }
               }

               /* if f(A)v = 0 for all generators v, we are done */
               if (i == n)
               {
                  nmod_mat_clear(mat);
                  nmod_poly_clear(poly);
                  break;
               }
            }

            nmod_mat_clear(mat);
            nmod_poly_clear(poly);
        }

        flint_free(P);
        flint_free(Q);
        fmpz_mat_clear(v2);
        fmpz_mat_clear(v1);
        fmpz_mat_clear(v3);
        fmpz_clear(m);
        _fmpz_vec_clear(rold, n);
    }

    return len;
}

void fmpz_mat_minpoly_modular(fmpz_poly_t cp, const fmpz_mat_t mat)
{
    slong len;

    fmpz_poly_fit_length(cp, mat->r + 1);

    len = _fmpz_mat_minpoly_modular(cp->coeffs, mat);

    _fmpz_poly_set_length(cp, len);
}