flint-sys 0.9.0

/*
    Copyright (C) 2023 Fredrik Johansson

    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 "fmpz.h"
#include "fmpz_factor.h"
#include "fmpz_vec.h"
#include "fmpz_mod.h"
#include "fmpz_mod_vec.h"
#include "fmpz_mod_mat.h"
#include "fmpz_mod_poly.h"
#include "fmpz_mod_poly_factor.h"
#include "gr.h"
#include "gr_vec.h"
#include "gr_poly.h"
#include "gr_mat.h"

typedef struct
{
    fmpz_mod_ctx_struct * ctx;
    truth_t is_prime;
    fmpz a;    /* when used as finite field with defining polynomial x - a */
}
_gr_fmpz_mod_ctx_struct;

#define FMPZ_MOD_CTX(ring_ctx) ((((_gr_fmpz_mod_ctx_struct *)(ring_ctx))->ctx))
#define FMPZ_MOD_IS_PRIME(ring_ctx) (((_gr_fmpz_mod_ctx_struct *)(ring_ctx))->is_prime)
#define FMPZ_MOD_CTX_A(ring_ctx) (&((((_gr_fmpz_mod_ctx_struct *)(ring_ctx))->a)))

static int
_gr_fmpz_mod_ctx_write(gr_stream_t out, gr_ctx_t ctx)
{
    gr_stream_write(out, "Integers mod ");
    gr_stream_write_fmpz(out, FMPZ_MOD_CTX(ctx)->n);
    gr_stream_write(out, " (fmpz)");
    return GR_SUCCESS;
}

static void
_gr_fmpz_mod_ctx_clear(gr_ctx_t ctx)
{
    fmpz_mod_ctx_clear(FMPZ_MOD_CTX(ctx));
    flint_free(FMPZ_MOD_CTX(ctx));
    fmpz_clear(FMPZ_MOD_CTX_A(ctx));
}

static int
_gr_fmpz_mod_ctx_set_is_field(gr_ctx_t ctx, truth_t is_field)
{
    FMPZ_MOD_IS_PRIME(ctx) = is_field;
    return GR_SUCCESS;
}

static truth_t
_gr_fmpz_mod_ctx_is_field(gr_ctx_t ctx)
{
/*
    if (FMPZ_MOD_IS_PRIME(ctx) == T_UNKNOWN)
        FMPZ_MOD_IS_PRIME(ctx) = fmpz_is_prime(FMPZ_MOD_CTX(ctx)->n) ? T_TRUE : T_FALSE;
*/

    return FMPZ_MOD_IS_PRIME(ctx);
}

static void
_gr_fmpz_mod_init(fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_init(x);
}

static void
_gr_fmpz_mod_clear(fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_clear(x);
}

static void
_gr_fmpz_mod_swap(fmpz_t x, fmpz_t y, const gr_ctx_t ctx)
{
    fmpz_t t;
    *t = *x;
    *x = *y;
    *y = *t;
}

static void
_gr_fmpz_mod_set_shallow(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    *res = *x;
}

static int
_gr_fmpz_mod_randtest(fmpz_t res, flint_rand_t state, const gr_ctx_t ctx)
{
    fmpz_mod_rand(res, state, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_write(gr_stream_t out, const fmpz_t x, const gr_ctx_t ctx)
{
    gr_stream_write_fmpz(out, x);
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_zero(fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_zero(x);
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_one(fmpz_t x, const gr_ctx_t ctx)
{
    if (fmpz_is_one(FMPZ_MOD_CTX(ctx)->n))
        fmpz_zero(x);
    else
        fmpz_one(x);
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_set_si(fmpz_t res, slong v, const gr_ctx_t ctx)
{
    fmpz_mod_set_si(res, v, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_set_ui(fmpz_t res, ulong v, const gr_ctx_t ctx)
{
    fmpz_mod_set_ui(res, v, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_get_fmpz(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_set(res, x);
    return GR_SUCCESS;
}


/* todo: public interface */
#define NMOD_CTX_REF(ring_ctx) (((nmod_t *)((ring_ctx))))
#define NMOD_CTX(ring_ctx) (*NMOD_CTX_REF(ring_ctx))

static int
_gr_fmpz_mod_set_other(fmpz_t res, gr_ptr v, gr_ctx_t v_ctx, const gr_ctx_t ctx)
{
    if (v_ctx->which_ring == GR_CTX_FMPZ_MOD)
    {
        if (!fmpz_equal(FMPZ_MOD_CTX(ctx)->n, FMPZ_MOD_CTX(v_ctx)->n))
            return GR_DOMAIN;

        fmpz_set(res, v);
        return GR_SUCCESS;
    }

    if (v_ctx->which_ring == GR_CTX_NMOD)
    {
        if (!fmpz_equal_ui(FMPZ_MOD_CTX(ctx)->n, NMOD_CTX(v_ctx).n))
            return GR_DOMAIN;

        fmpz_set_ui(res, * (ulong *) v);
        return GR_SUCCESS;
    }

    return GR_UNABLE;
}

static int
_gr_fmpz_mod_set_fmpz(fmpz_t res, const fmpz_t v, const gr_ctx_t ctx)
{
    fmpz_mod_set_fmpz(res, v, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static truth_t
_gr_fmpz_mod_is_zero(const fmpz_t x, const gr_ctx_t ctx)
{
    return fmpz_is_zero(x) ? T_TRUE : T_FALSE;
}

static truth_t
_gr_fmpz_mod_is_one(const fmpz_t x, const gr_ctx_t ctx)
{
    return fmpz_mod_is_one(x, FMPZ_MOD_CTX(ctx)) ? T_TRUE : T_FALSE;
}

static truth_t
_gr_fmpz_mod_is_neg_one(const fmpz_t x, const gr_ctx_t ctx)
{
    truth_t res;
    fmpz_t t;
    fmpz_init(t);
    fmpz_mod_set_si(t, -1, FMPZ_MOD_CTX(ctx));
    res = fmpz_equal(t, x) ? T_TRUE : T_FALSE;
    fmpz_clear(t);
    return res;
}

static truth_t
_gr_fmpz_mod_equal(const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    return fmpz_equal(x, y) ? T_TRUE : T_FALSE;
}

static int
_gr_fmpz_mod_set(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_set(res, x);
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_neg(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_mod_neg(res, x, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_add(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    fmpz_mod_add(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_add_si(fmpz_t res, const fmpz_t x, slong y, const gr_ctx_t ctx)
{
    fmpz_mod_add_si(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_add_ui(fmpz_t res, const fmpz_t x, ulong y, const gr_ctx_t ctx)
{
    fmpz_mod_add_ui(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_sub(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    fmpz_mod_sub(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_sub_si(fmpz_t res, const fmpz_t x, slong y, const gr_ctx_t ctx)
{
    fmpz_mod_sub_si(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_sub_ui(fmpz_t res, const fmpz_t x, ulong y, const gr_ctx_t ctx)
{
    fmpz_mod_sub_ui(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_mul(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
#if 1
    fmpz_mod_mul(res, x, y, FMPZ_MOD_CTX(ctx));
#else
    fmpz_mul(res, x, y);
    fmpz_mod(res, res, FMPZ_MOD_CTX(ctx)->n);
#endif
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_mul_si(fmpz_t res, const fmpz_t x, slong y, const gr_ctx_t ctx)
{
    fmpz_mod_mul_si(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_mul_ui(fmpz_t res, const fmpz_t x, ulong y, const gr_ctx_t ctx)
{
    fmpz_mod_mul_ui(res, x, y, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_addmul(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    fmpz_t t;
    fmpz_init(t);
    fmpz_mul(t, x, y);
    fmpz_add(t, t, res);
    fmpz_mod_set_fmpz(res, t, FMPZ_MOD_CTX(ctx));
    fmpz_clear(t);
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_submul(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    fmpz_t t;
    fmpz_init(t);
    fmpz_mul(t, x, y);
    fmpz_sub(t, res, t);
    fmpz_mod_set_fmpz(res, t, FMPZ_MOD_CTX(ctx));
    fmpz_clear(t);
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_mul_two(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_mod_add(res, x, x, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_sqr(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    fmpz_mod_mul(res, x, x, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_inv(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    /* todo: also check for -1 when fast? */
    if (fmpz_is_one(x))
    {
        fmpz_one(res);
        return GR_SUCCESS;
    }
    else if (fmpz_is_zero(x))
    {
        fmpz_zero(res);
        return fmpz_is_one(FMPZ_MOD_CTX(ctx)->n) ? GR_SUCCESS : GR_DOMAIN;
    }
    else
    {
        int status;

        fmpz_t d;
        fmpz_init(d);
        fmpz_gcdinv(d, res, x, FMPZ_MOD_CTX(ctx)->n);

        if (fmpz_is_one(d))
            status = GR_SUCCESS;
        else
            status = GR_DOMAIN;

        fmpz_clear(d);
        return status;
    }
}

static int
_gr_fmpz_mod_div(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    int status;
    fmpz_t t;

    fmpz_init(t);
    status = _gr_fmpz_mod_inv(t, y, ctx);
    if (status == GR_SUCCESS)
        fmpz_mod_mul(res, x, t, FMPZ_MOD_CTX(ctx));
    else
        fmpz_zero(res);
    fmpz_clear(t);

    return status;
}

static int
_gr_fmpz_mod_div_nonunique(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    int status;

#if 1
    status = fmpz_mod_divides(res, x, y, FMPZ_MOD_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
#else
    if (FMPZ_MOD_IS_PRIME(ctx) != T_TRUE)
    {
        status = fmpz_mod_divides(res, x, y, FMPZ_MOD_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
    }
    else
    {
        fmpz_t t;
        fmpz_init(t);
        status = _gr_fmpz_mod_inv(t, y, ctx);
        if (status == GR_SUCCESS)
            fmpz_mod_mul(res, x, t, FMPZ_MOD_CTX(ctx));
        else
            fmpz_zero(res);
        fmpz_clear(t);
    }
#endif

    return status;
}

static truth_t
_gr_fmpz_mod_divides(const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
    truth_t res;
    fmpz_t t;
    fmpz_init(t);
    res = fmpz_mod_divides(t, y, x, FMPZ_MOD_CTX(ctx)) ? T_TRUE : T_FALSE;
    fmpz_clear(t);
    return res;
}

static truth_t
_gr_fmpz_mod_is_invertible(const fmpz_t x, const gr_ctx_t ctx)
{
    return fmpz_mod_is_invertible(x, FMPZ_MOD_CTX(ctx)) ? T_TRUE : T_FALSE;
}

static int
_gr_fmpz_mod_pow_ui(fmpz_t res, const fmpz_t x, ulong exp, const gr_ctx_t ctx)
{
    fmpz_mod_pow_ui(res, x, exp, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_pow_fmpz(fmpz_t res, const fmpz_t x, const fmpz_t exp, const gr_ctx_t ctx)
{
    return fmpz_mod_pow_fmpz(res, x, exp, FMPZ_MOD_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
}

static int
_gr_fmpz_mod_sqrt(fmpz_t res, const fmpz_t x, const gr_ctx_t ctx)
{
    if (fmpz_is_zero(x) || fmpz_is_one(x))
    {
        fmpz_set(res, x);
        return GR_SUCCESS;
    }
    else if (FMPZ_MOD_IS_PRIME(ctx) == T_TRUE)
    {
        return fmpz_sqrtmod(res, x, FMPZ_MOD_CTX(ctx)->n) ? GR_SUCCESS : GR_DOMAIN;
    }
    else
    {
        /* todo: implement the general case */
        return GR_UNABLE;
    }
}

static truth_t
_gr_fmpz_mod_is_square(const fmpz_t x, const gr_ctx_t ctx)
{
    if (fmpz_is_zero(x) || fmpz_is_one(x))
    {
        return T_TRUE;
    }
    else if (FMPZ_MOD_IS_PRIME(ctx) == T_TRUE)
    {
        fmpz_t t;
        truth_t ans;
        fmpz_init(t);
        ans = fmpz_sqrtmod(t, x, FMPZ_MOD_CTX(ctx)->n) ? T_TRUE : T_FALSE;
        fmpz_clear(t);
        return ans;
    }
    else
    {
        /* todo: implement the general case */
        return T_UNKNOWN;
    }
}

/* todo: len 1 */
static int
_gr_fmpz_mod_vec_dot(fmpz_t res, const fmpz_t initial, int subtract, const fmpz * vec1, const fmpz * vec2, slong len, gr_ctx_t ctx)
{
    if (len <= 0)
    {
        if (initial == NULL)
            fmpz_zero(res);
        else
            fmpz_set(res, initial);
        return GR_SUCCESS;
    }

    _fmpz_vec_dot_general(res, initial, subtract, vec1, vec2, 0, len);
    fmpz_mod_set_fmpz(res, res, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

/* todo: len 1 */
static int
_gr_fmpz_mod_vec_dot_rev(fmpz_t res, const fmpz_t initial, int subtract, const fmpz * vec1, const fmpz * vec2, slong len, gr_ctx_t ctx)
{
    if (len <= 0)
    {
        if (initial == NULL)
            fmpz_zero(res);
        else
            fmpz_set(res, initial);
        return GR_SUCCESS;
    }

    _fmpz_vec_dot_general(res, initial, subtract, vec1, vec2, 1, len);
    fmpz_mod_set_fmpz(res, res, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_vec_mul_scalar(fmpz * res, const fmpz * vec, slong len, const fmpz_t c, gr_ctx_t ctx)
{
    _fmpz_mod_vec_scalar_mul_fmpz_mod(res, vec, len, c, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_scalar_mul_vec(fmpz * res, fmpz_t c, const fmpz * vec, slong len, gr_ctx_t ctx)
{
    _fmpz_mod_vec_scalar_mul_fmpz_mod(res, vec, len, c, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_vec_addmul_scalar(fmpz * res, const fmpz * vec, slong len, const fmpz_t c, gr_ctx_t ctx)
{
    _fmpz_mod_vec_scalar_addmul_fmpz_mod(res, vec, len, c, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}


static int
_gr_fmpz_mod_poly_mullow(fmpz * res,
    const fmpz * poly1, slong len1,
    const fmpz * poly2, slong len2, slong n, gr_ctx_t ctx)
{
    if (len1 >= len2)
        _fmpz_mod_poly_mullow(res, poly1, len1, poly2, len2, n, FMPZ_MOD_CTX(ctx));
    else
        _fmpz_mod_poly_mullow(res, poly2, len2, poly1, len1, n, FMPZ_MOD_CTX(ctx));

    return GR_SUCCESS;
}

static int
_gr_fmpz_mod_poly_mulmid(fmpz * res,
    const fmpz * poly1, slong len1,
    const fmpz * poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx)
{
    if (len1 >= len2)
        _fmpz_mod_poly_mulmid(res, poly1, len1, poly2, len2, nlo, nhi, FMPZ_MOD_CTX(ctx));
    else
        _fmpz_mod_poly_mulmid(res, poly2, len2, poly1, len1, nlo, nhi, FMPZ_MOD_CTX(ctx));

    return GR_SUCCESS;
}


/* fixme: duplicate _fmpz_mod_poly methods for error handling */
static int
_gr_fmpz_mod_poly_divrem(fmpz * Q, fmpz * R, const fmpz * A, slong lenA,
                                  const fmpz * B, slong lenB, gr_ctx_t ctx)
{
    if (lenB <= 30 || lenA - lenB <= 5)
    {
        fmpz_t invB;
        int status;

        fmpz_init(invB);
        status = _gr_fmpz_mod_inv(invB, B + lenB - 1, ctx);
        if (status == GR_SUCCESS)
            _fmpz_mod_poly_divrem_basecase(Q, R, A, lenA, B, lenB, invB, FMPZ_MOD_CTX(ctx));

        fmpz_clear(invB);
        return status;
    }
    else
    {
        return _gr_poly_divrem_newton(Q, R, A, lenA, B, lenB, ctx);
    }
}

#define TUNE_TAB_SIZE 23

static const int tuning_bit_steps[TUNE_TAB_SIZE] = { 32, 45, 64, 91, 128, 181, 256, 362, 512, 724, 1024, 1448, 2048, 2896, 4096, 5793, 8192, 11585, 16384, 23170, 32768, 46341, 65536};
static const short inv_series_cutoff_tab[TUNE_TAB_SIZE] = {21, 14, 40, 39, 48, 60, 89, 72, 72, 54, 48, 39, 32, 24, 24, 20, 17, 18, 16, 15, 13, 12, 14, };
static const short div_series_cutoff_tab[TUNE_TAB_SIZE] = {23, 21, 52, 50, 66, 101, 106, 97, 106, 72, 60, 50, 44, 35, 38, 30, 26, 22, 20, 18, 16, 14, 22, };

static const slong find_cutoff(const short * tab, slong b)
{
    slong i;

    i = 0;
    while (i + 1 < TUNE_TAB_SIZE && tuning_bit_steps[i + 1] <= b)
        i++;

    return tab[i];
}

static int
_gr_fmpz_mod_poly_inv_series(fmpz * Q, const fmpz * B, slong lenB, slong len, gr_ctx_t ctx)
{
    slong cutoff, bits;
    lenB = FLINT_MIN(len, lenB);

    if (lenB <= 20)
        return _gr_poly_inv_series_basecase(Q, B, lenB, len, ctx);

    bits = fmpz_bits(fmpz_mod_ctx_modulus(FMPZ_MOD_CTX(ctx)));
    cutoff  = find_cutoff(inv_series_cutoff_tab, bits);

    if (lenB <= cutoff)
        return _gr_poly_inv_series_basecase(Q, B, lenB, len, ctx);
    else
        return _gr_poly_inv_series_newton(Q, B, lenB, len, cutoff, ctx);
}

/* todo: the fmpz_mod_poly module has better basecase code */
static int
_gr_fmpz_mod_poly_div_series(fmpz * Q, const fmpz * A, slong lenA, const fmpz * B, slong lenB, slong len, gr_ctx_t ctx)
{
    slong cutoff, bits;

    lenA = FLINT_MIN(len, lenA);
    lenB = FLINT_MIN(len, lenB);

    if (lenB <= 20)
        return _gr_poly_div_series_basecase(Q, A, lenA, B, lenB, len, ctx);

    bits = fmpz_bits(fmpz_mod_ctx_modulus(FMPZ_MOD_CTX(ctx)));
    cutoff  = find_cutoff(div_series_cutoff_tab, bits);

    if (lenB <= cutoff)
        return _gr_poly_div_series_basecase(Q, A, lenA, B, lenB, len, ctx);
    else
        return _gr_poly_div_series_newton(Q, A, lenA, B, lenB, len, cutoff, ctx);
}

static int _gr_fmpz_mod_poly_gcd(nn_ptr G, slong * lenG, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx)
{
    if (FLINT_MIN(lenA, lenB) < FMPZ_MOD_POLY_GCD_CUTOFF)
        return _gr_poly_gcd_euclidean(G, lenG, A, lenA, B, lenB, ctx);
    else
        return _gr_poly_gcd_hgcd(G, lenG, A, lenA, B, lenB, FMPZ_MOD_POLY_HGCD_CUTOFF, FMPZ_MOD_POLY_GCD_CUTOFF, ctx);
}


/* todo: also need the _other version ... ? */
/* todo: implement generically */

static int
_gr_fmpz_mod_roots_gr_poly(gr_vec_t roots, gr_vec_t mult, const fmpz_mod_poly_t poly, int flags, gr_ctx_t ctx)
{
    if (poly->length == 0)
        return GR_DOMAIN;

    {
        gr_ctx_t ZZ;
        fmpz_mod_poly_factor_t fac;
        slong i, num;
        int status = GR_SUCCESS;

        gr_ctx_init_fmpz(ZZ);

        fmpz_mod_poly_factor_init(fac, FMPZ_MOD_CTX(ctx));

        if (gr_ctx_is_field(ctx) == T_TRUE)
        {
            fmpz_mod_poly_roots(fac, poly, 1, FMPZ_MOD_CTX(ctx));
        }
        else
        {
            fmpz_factor_t nfac;
            fmpz_factor_init(nfac);
            fmpz_factor(nfac, FMPZ_MOD_CTX(ctx)->n);

            num = 0;
            for (i = 0; i < nfac->num; i++)
                num += nfac->exp[i];

            if (num > 20)
            {
                status = GR_UNABLE;
            }
            else
            {
                if (!fmpz_mod_poly_roots_factored_with_length_limit(fac, poly, 1, 1000000, nfac, FMPZ_MOD_CTX(ctx)))
                    status = GR_UNABLE;
            }

            fmpz_factor_clear(nfac);
        }

        if (status == GR_SUCCESS)
        {
            num = fac->num;

            gr_vec_set_length(roots, num, ctx);
            gr_vec_set_length(mult, num, ZZ);

            for (i = 0; i < num; i++)
            {
                fmpz_mod_neg(gr_vec_entry_ptr(roots, i, ctx), fac->poly[i].coeffs, FMPZ_MOD_CTX(ctx));

                /* work around flint bug: factors can be non-monic */
                if (!fmpz_mod_is_one(fac->poly[i].coeffs + 1, FMPZ_MOD_CTX(ctx)))
                    status |= _gr_fmpz_mod_div(gr_vec_entry_ptr(roots, i, ctx), gr_vec_entry_ptr(roots, i, ctx), fac->poly[i].coeffs + 1, ctx);

                fmpz_set_ui(((fmpz *) mult->entries) + i, fac->exp[i]);
            }
        }

        fmpz_mod_poly_factor_clear(fac, FMPZ_MOD_CTX(ctx));
        gr_ctx_clear(ZZ);

        return status;
    }
}

static int
_gr_fmpz_mod_mat_mul(fmpz_mod_mat_t res, const fmpz_mod_mat_t x, const fmpz_mod_mat_t y, gr_ctx_t ctx)
{
    fmpz_mat_mul(res, x, y);
    _fmpz_mod_mat_reduce(res, FMPZ_MOD_CTX(ctx));
    return GR_SUCCESS;
}

/* todo: tune cutoff for different bit sizes */
/* also tune cutoff for triangular solving */
static int
_gr_fmpz_mod_mat_lu(slong * rank, slong * P, fmpz_mod_mat_t LU, const fmpz_mod_mat_t A, int rank_check, gr_ctx_t ctx)
{
    slong cutoff = 8;

    if (A->r < cutoff || A->c < cutoff)
        return gr_mat_lu_classical(rank, P, (gr_mat_struct *) LU, (const gr_mat_struct *) A, rank_check, ctx);
    else
        return gr_mat_lu_recursive(rank, P, (gr_mat_struct *) LU, (const gr_mat_struct *) A, rank_check, ctx);
}

static int
_gr_fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, gr_ctx_t ctx)
{
    slong n = mat->r;

    if (n <= 4)
        return gr_mat_det_cofactor(res, (const gr_mat_struct *) mat, ctx);

    if (gr_mat_det_lu(res, (const gr_mat_struct *) mat, ctx) == GR_SUCCESS)
        return GR_SUCCESS;

    /* Fall back on division-free algorithm if we encountered an impossible inverse */
    /* Could try something else here: faddeev_bsgs (O(n^3.5)) or Howell form. */
    return gr_mat_det_berkowitz(res, (const gr_mat_struct *) mat, ctx);
}

static int
_gr_fmpz_mod_mat_charpoly(fmpz * res, const fmpz_mod_mat_t mat, gr_ctx_t ctx)
{
    slong n = mat->r;

    if (n >= 16 && _gr_mat_charpoly_danilevsky(res, (const gr_mat_struct *) mat, ctx) == GR_SUCCESS)
        return GR_SUCCESS;

    return _gr_mat_charpoly_berkowitz(res, (const gr_mat_struct *) mat, ctx);
}

static int
_gr_fmpz_mod_mat_reduce_row(slong * column, fmpz_mod_mat_t mat, slong * P, slong * L, slong n, gr_ctx_t ctx)
{
    return fmpz_mod_mat_reduce_row(column, mat, P, L, n, FMPZ_MOD_CTX(ctx));
}


int _fmpz_mod_methods_initialized = 0;

gr_static_method_table _fmpz_mod_methods;

gr_method_tab_input _fmpz_mod_methods_input[] =
{
    {GR_METHOD_CTX_WRITE,       (gr_funcptr) _gr_fmpz_mod_ctx_write},
    {GR_METHOD_CTX_CLEAR,       (gr_funcptr) _gr_fmpz_mod_ctx_clear},
    {GR_METHOD_CTX_IS_RING,     (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_COMMUTATIVE_RING, (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_INTEGRAL_DOMAIN,  (gr_funcptr) _gr_fmpz_mod_ctx_is_field},
    {GR_METHOD_CTX_IS_FIELD,            (gr_funcptr) _gr_fmpz_mod_ctx_is_field},
    {GR_METHOD_CTX_IS_UNIQUE_FACTORIZATION_DOMAIN,
                                (gr_funcptr) _gr_fmpz_mod_ctx_is_field},
    {GR_METHOD_CTX_IS_FINITE,
                                (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_FINITE_CHARACTERISTIC,
                                (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_EXACT,    (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_CANONICAL,
                                (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_SET_IS_FIELD,(gr_funcptr) _gr_fmpz_mod_ctx_set_is_field},
    {GR_METHOD_INIT,            (gr_funcptr) _gr_fmpz_mod_init},
    {GR_METHOD_CLEAR,           (gr_funcptr) _gr_fmpz_mod_clear},
    {GR_METHOD_SWAP,            (gr_funcptr) _gr_fmpz_mod_swap},
    {GR_METHOD_SET_SHALLOW,     (gr_funcptr) _gr_fmpz_mod_set_shallow},
    {GR_METHOD_RANDTEST,        (gr_funcptr) _gr_fmpz_mod_randtest},
    {GR_METHOD_WRITE,           (gr_funcptr) _gr_fmpz_mod_write},
    {GR_METHOD_ZERO,            (gr_funcptr) _gr_fmpz_mod_zero},
    {GR_METHOD_ONE,             (gr_funcptr) _gr_fmpz_mod_one},
    {GR_METHOD_IS_ZERO,         (gr_funcptr) _gr_fmpz_mod_is_zero},
    {GR_METHOD_IS_ONE,          (gr_funcptr) _gr_fmpz_mod_is_one},
    {GR_METHOD_IS_NEG_ONE,      (gr_funcptr) _gr_fmpz_mod_is_neg_one},
    {GR_METHOD_EQUAL,           (gr_funcptr) _gr_fmpz_mod_equal},
    {GR_METHOD_SET,             (gr_funcptr) _gr_fmpz_mod_set},
    {GR_METHOD_SET_SI,          (gr_funcptr) _gr_fmpz_mod_set_si},
    {GR_METHOD_SET_UI,          (gr_funcptr) _gr_fmpz_mod_set_ui},
    {GR_METHOD_SET_FMPZ,        (gr_funcptr) _gr_fmpz_mod_set_fmpz},
    {GR_METHOD_SET_OTHER,       (gr_funcptr) _gr_fmpz_mod_set_other},
    {GR_METHOD_GET_FMPZ,        (gr_funcptr) _gr_fmpz_mod_get_fmpz},
    {GR_METHOD_NEG,             (gr_funcptr) _gr_fmpz_mod_neg},
    {GR_METHOD_ADD,             (gr_funcptr) _gr_fmpz_mod_add},
    {GR_METHOD_ADD_UI,          (gr_funcptr) _gr_fmpz_mod_add_ui},
    {GR_METHOD_ADD_SI,          (gr_funcptr) _gr_fmpz_mod_add_si},
    {GR_METHOD_SUB,             (gr_funcptr) _gr_fmpz_mod_sub},
    {GR_METHOD_SUB_UI,          (gr_funcptr) _gr_fmpz_mod_sub_ui},
    {GR_METHOD_SUB_SI,          (gr_funcptr) _gr_fmpz_mod_sub_si},
    {GR_METHOD_MUL,             (gr_funcptr) _gr_fmpz_mod_mul},
    {GR_METHOD_MUL_SI,          (gr_funcptr) _gr_fmpz_mod_mul_si},
    {GR_METHOD_MUL_UI,          (gr_funcptr) _gr_fmpz_mod_mul_ui},
    {GR_METHOD_ADDMUL,          (gr_funcptr) _gr_fmpz_mod_addmul},
    {GR_METHOD_SUBMUL,          (gr_funcptr) _gr_fmpz_mod_submul},
    {GR_METHOD_MUL_TWO,         (gr_funcptr) _gr_fmpz_mod_mul_two},
    {GR_METHOD_SQR,             (gr_funcptr) _gr_fmpz_mod_sqr},
    {GR_METHOD_DIV,             (gr_funcptr) _gr_fmpz_mod_div},
    {GR_METHOD_DIV_NONUNIQUE,   (gr_funcptr) _gr_fmpz_mod_div_nonunique},
    {GR_METHOD_DIVIDES,         (gr_funcptr) _gr_fmpz_mod_divides},
    {GR_METHOD_IS_INVERTIBLE,   (gr_funcptr) _gr_fmpz_mod_is_invertible},
    {GR_METHOD_INV,             (gr_funcptr) _gr_fmpz_mod_inv},
    {GR_METHOD_POW_UI,          (gr_funcptr) _gr_fmpz_mod_pow_ui},
    {GR_METHOD_POW_FMPZ,        (gr_funcptr) _gr_fmpz_mod_pow_fmpz},
    {GR_METHOD_SQRT,            (gr_funcptr) _gr_fmpz_mod_sqrt},
    {GR_METHOD_IS_SQUARE,       (gr_funcptr) _gr_fmpz_mod_is_square},

/*
    {GR_METHOD_VEC_INIT,        (gr_funcptr) _gr_mpn_mod_vec_zero},
    {GR_METHOD_VEC_CLEAR,       (gr_funcptr) _gr_mpn_mod_vec_clear},
    {GR_METHOD_VEC_SET,         (gr_funcptr) _gr_mpn_mod_vec_set},
    {GR_METHOD_VEC_SWAP,        (gr_funcptr) _gr_mpn_mod_vec_swap},
    {GR_METHOD_VEC_ZERO,        (gr_funcptr) _gr_mpn_mod_vec_zero},
    {GR_METHOD_VEC_NEG,         (gr_funcptr) _gr_mpn_mod_vec_neg},
    {GR_METHOD_VEC_ADD,         (gr_funcptr) _gr_mpn_mod_vec_add},
    {GR_METHOD_VEC_SUB,         (gr_funcptr) _gr_mpn_mod_vec_sub},
    {GR_METHOD_VEC_MUL,         (gr_funcptr) _gr_mpn_mod_vec_mul},
*/
    {GR_METHOD_VEC_MUL_SCALAR,  (gr_funcptr) _gr_fmpz_mod_vec_mul_scalar},
    {GR_METHOD_SCALAR_MUL_VEC,  (gr_funcptr) _gr_fmpz_mod_scalar_mul_vec},
    {GR_METHOD_VEC_ADDMUL_SCALAR,    (gr_funcptr) _gr_fmpz_mod_vec_addmul_scalar},

    {GR_METHOD_VEC_DOT,         (gr_funcptr) _gr_fmpz_mod_vec_dot},
    {GR_METHOD_VEC_DOT_REV,     (gr_funcptr) _gr_fmpz_mod_vec_dot_rev},
    {GR_METHOD_POLY_MULLOW,     (gr_funcptr) _gr_fmpz_mod_poly_mullow},
    {GR_METHOD_POLY_MULMID,     (gr_funcptr) _gr_fmpz_mod_poly_mulmid},
    {GR_METHOD_POLY_INV_SERIES, (gr_funcptr) _gr_fmpz_mod_poly_inv_series},
    {GR_METHOD_POLY_DIV_SERIES, (gr_funcptr) _gr_fmpz_mod_poly_div_series},
    {GR_METHOD_POLY_DIVREM,     (gr_funcptr) _gr_fmpz_mod_poly_divrem},
    {GR_METHOD_POLY_GCD,        (gr_funcptr) _gr_fmpz_mod_poly_gcd},
    {GR_METHOD_POLY_ROOTS,      (gr_funcptr) _gr_fmpz_mod_roots_gr_poly},
    {GR_METHOD_MAT_MUL,         (gr_funcptr) _gr_fmpz_mod_mat_mul},
    {GR_METHOD_MAT_LU,          (gr_funcptr) _gr_fmpz_mod_mat_lu},
    {GR_METHOD_MAT_DET,         (gr_funcptr) _gr_fmpz_mod_mat_det},
    {GR_METHOD_MAT_CHARPOLY,    (gr_funcptr) _gr_fmpz_mod_mat_charpoly},
    {GR_METHOD_MAT_REDUCE_ROW,  (gr_funcptr) _gr_fmpz_mod_mat_reduce_row},
    {0,                         (gr_funcptr) NULL},
};

void
gr_ctx_init_fmpz_mod(gr_ctx_t ctx, const fmpz_t n)
{
    ctx->which_ring = GR_CTX_FMPZ_MOD;
    ctx->sizeof_elem = sizeof(fmpz);

    FMPZ_MOD_CTX(ctx) = flint_malloc(sizeof(fmpz_mod_ctx_struct));
    fmpz_mod_ctx_init(FMPZ_MOD_CTX(ctx), n);
    FMPZ_MOD_IS_PRIME(ctx) = T_UNKNOWN;
    fmpz_init(FMPZ_MOD_CTX_A(ctx));

    ctx->size_limit = WORD_MAX;

    ctx->methods = _fmpz_mod_methods;

    if (!_fmpz_mod_methods_initialized)
    {
        gr_method_tab_init(_fmpz_mod_methods, _fmpz_mod_methods_input);
        _fmpz_mod_methods_initialized = 1;
    }
}

void
_gr_ctx_init_fmpz_mod_from_ref(gr_ctx_t ctx, const void * fctx)
{
    ctx->which_ring = GR_CTX_FMPZ_MOD;
    ctx->sizeof_elem = sizeof(fmpz);

    FMPZ_MOD_CTX(ctx) = (fmpz_mod_ctx_struct *) fctx;
    FMPZ_MOD_IS_PRIME(ctx) = T_UNKNOWN;
    fmpz_init(FMPZ_MOD_CTX_A(ctx));

    ctx->size_limit = WORD_MAX;

    ctx->methods = _fmpz_mod_methods;

    if (!_fmpz_mod_methods_initialized)
    {
        gr_method_tab_init(_fmpz_mod_methods, _fmpz_mod_methods_input);
        _fmpz_mod_methods_initialized = 1;
    }
}