flint-sys 0.9.0

/*
    Copyright (C) 2007 David Howden
    Copyright (C) 2007-2010, 2020, 2022 William Hart
    Copyright (C) 2008 Richard Howell-Peak
    Copyright (C) 2011 Fredrik Johansson
    Copyright (C) 2012 Lina Kulakova
    Copyright (C) 2013 Martin Lee

    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 "nmod_poly_factor.h"

/* Return 1 if poly has a trivial factor, otherwise return 0. For small p,
   this quickly filters out many candidates when testing random polynomials
   for irreducibility. */
static int nmod_poly_is_reducible_trial_div(const nmod_poly_t poly)
{
    ulong x, p = poly->mod.n;

    if (poly->mod.n > FLINT_MAX(200, 2 * poly->length))
        return 0;

    /* Try all linear factors. To do: check if it is worthwhile to eliminate
       other low-degree factors too. */
    /* To do: use multipoint evaluation when p is large enough. */
    /* To do: use sparse algorithm when poly is sparse. */
    for (x = 1; x < p; x++)
        if (nmod_poly_evaluate_nmod(poly, x) == 0)
            return 1;

    return 0;
}

int nmod_poly_is_irreducible_ddf(const nmod_poly_t poly)
{

    nmod_poly_t f, v, vinv, tmp;
    nmod_poly_struct * h, * H, * I;
    nmod_mat_t HH;
    slong i, j, l, m, n, d;
    double beta;
    int result = 1;
    n = nmod_poly_degree(poly);

    if (n < 2)
        return 1;

    if (!nmod_poly_is_squarefree(poly))
        return 0;

    beta = 0.5 * (1. - (log(2)/log(n)));
    l = ceil(pow (n, beta));
    m = ceil(0.5*n/l);

    /* initialization */
    nmod_poly_init_mod(f, poly->mod);
    nmod_poly_init_mod(v, poly->mod);
    nmod_poly_init_mod(vinv, poly->mod);
    nmod_poly_init_mod(tmp, poly->mod);

    h =  flint_malloc((2 * m + l + 1) * sizeof(nmod_poly_struct));
    H = h + (l + 1);
    I = H + m;

    for (i = 0; i < 2*m + l + 1; i++)
        nmod_poly_init_mod(h + i, poly->mod);

    nmod_poly_make_monic(v, poly);

    nmod_poly_reverse(vinv, v, v->length);
    nmod_poly_inv_series(vinv, vinv, v->length);

    /* compute baby steps: h[i] = x^{p^i}mod v */
    nmod_poly_set_coeff_ui(h + 0, 1, 1);
    nmod_poly_powmod_x_ui_preinv(h + 1, poly->mod.n, v, vinv);

    if (FLINT_BIT_COUNT(poly->mod.n) > ((n_sqrt(v->length - 1) + 1)*3)/4)
    {
        for (i = 1; i < (slong) FLINT_BIT_COUNT(l); i++)
            nmod_poly_compose_mod_brent_kung_vec_preinv(h + 1 +
                            (1 << (i - 1)), h + 1, (1 << (i - 1)),
                            (1 << (i - 1)), h + (1 << (i - 1)), v, vinv);

        nmod_poly_compose_mod_brent_kung_vec_preinv(h + 1 + (1 << (i - 1)),
                            h + 1, (1 << (i - 1)), l - (1 << (i - 1)),
						    h + (1 << (i - 1)), v, vinv);
    }
    else
    {
        for (i = 2; i < l + 1; i++)
        {
            nmod_poly_init_mod(h + i, poly->mod);

            nmod_poly_powmod_ui_binexp_preinv(h + i, h + i - 1, poly->mod.n,
                                              v, vinv);
        }
    }

    /* compute coarse distinct-degree factorisation */
    nmod_poly_set(H + 0, h + l);
    nmod_mat_init(HH, n_sqrt(v->length - 1) + 1, v->length - 1, poly->mod.n);
    nmod_poly_precompute_matrix(HH, H + 0, v, vinv);

    d = 1;
    for (j = 0; j < m; j++)
    {
        /* compute giant steps: H[j] = x^{p^(lj)}mod s */
        if (j > 0)
            nmod_poly_compose_mod_brent_kung_precomp_preinv(H + j, H + j - 1, HH,
                                                            v, vinv);
        /* compute interval polynomials */
        nmod_poly_set_coeff_ui(I + j, 0, 1);

        for (i = l - 1; i >= 0 && 2*d <= v->length - 1; i--, d++)
        {
            nmod_poly_rem(tmp, h + i, v);
            nmod_poly_sub(tmp, H + j, tmp);
            nmod_poly_mulmod_preinv (I + j, tmp, I + j, v, vinv);
        }

        /* compute F_j=f^{[j*l+1]} * ... * f^{[j*l+l]} */
        /* F_j is stored on the place of I_j */
        nmod_poly_gcd(I + j, v, I + j);

        if (I[j].length > 1)
        {
            result = 0;
            break;
        }
    }

    nmod_poly_clear(f);
    nmod_poly_clear(v);
    nmod_poly_clear(vinv);
    nmod_poly_clear(tmp);

    nmod_mat_clear (HH);

    for (i = 0; i < l + 1; i++)
        nmod_poly_clear(h + i);

    for (i = 0; i < m; i++)
    {
        nmod_poly_clear(H + i);
        nmod_poly_clear(I + i);
    }

    flint_free (h);

    return result;
}

int
nmod_poly_is_irreducible(const nmod_poly_t f)
{
    if (nmod_poly_length(f) <= 2)
        return 1;

    if (nmod_poly_is_reducible_trial_div(f))
        return 0;

    return nmod_poly_is_irreducible_ddf(f);
}

static void
nmod_poly_powpowmod(nmod_poly_t res, const nmod_poly_t pol,
                                    ulong exp, ulong exp2, const nmod_poly_t f)
{
    nmod_poly_t pow;
    ulong i;

    nmod_poly_init_mod(pow, f->mod);

    nmod_poly_powmod_ui_binexp(pow, pol, exp, f);

    nmod_poly_set(res, pow);

    if (!nmod_poly_equal(pow, pol))
        for (i = 1; i < exp2; i++)
            nmod_poly_powmod_ui_binexp(res, res, exp, f);

    nmod_poly_clear(pow);
}

int
nmod_poly_is_irreducible_rabin(const nmod_poly_t f)
{
    if (nmod_poly_length(f) > 2)
    {
        const ulong p = nmod_poly_modulus(f);
        const slong n     = nmod_poly_degree(f);
        nmod_poly_t a, x, x_p;

        nmod_poly_init(a, p);
        nmod_poly_init(x, p);
        nmod_poly_init(x_p, p);

	    nmod_poly_set_coeff_ui(x, 1, 1);

        /* Compute x^q mod f */
        nmod_poly_powpowmod(x_p, x, p, n, f);

	    if (!nmod_poly_is_zero(x_p))
            nmod_poly_make_monic(x_p, x_p);

        /* Now do the irreducibility test */
        if (!nmod_poly_equal(x_p, x))
        {
            nmod_poly_clear(a);
            nmod_poly_clear(x);
            nmod_poly_clear(x_p);

	    return 0;
        } else
        {
            n_factor_t factors;
            slong i;

            n_factor_init(&factors);

	        n_factor(&factors, n, 1);

            for (i = 0; i < factors.num; i++)
            {
                nmod_poly_powpowmod(a, x, p, n/factors.p[i], f);
                nmod_poly_sub(a, a, x);

                if (!nmod_poly_is_zero(a))
                    nmod_poly_make_monic(a, a);

                nmod_poly_gcd(a, a, f);

                if (a->length != 1)
                {
                    nmod_poly_clear(a);
                    nmod_poly_clear(x);
                    nmod_poly_clear(x_p);

		    return 0;
                }
            }
        }

        nmod_poly_clear(a);
        nmod_poly_clear(x);
        nmod_poly_clear(x_p);
    }

    return 1;
}