flint-sys 0.9.0

Bindings to the FLINT C library
Documentation
/*
    Copyright (C) 2024 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 "mpn_mod.h"
#include "gr_poly.h"

/* generated by gr_poly/tune/cutoffs */

static const short inv_series_cutoff_tab[] = {200, 220, 220, 279, 127, 139, 166, 166, 166, 200, 200, 210, 166, 166, 166, 191, 166, 166, 200, 152, 139, 139, 166, 145, 139, 145, 145, 133, 106, 133, 133, 116, 116, 133, 133, 106, 106, 106, 133, 97, 78, 97, 133, 85, 72, 85, 101, 78, 72, 72, 78, 72, 72, 72, 72, 66, 66, 66, 66, 66, };
static const short div_series_cutoff_tab[] = {231, 306, 321, 370, 166, 182, 220, 210, 231, 242, 254, 254, 210, 220, 220, 210, 200, 191, 210, 210, 166, 166, 182, 191, 152, 145, 159, 166, 133, 133, 152, 145, 116, 133, 133, 127, 121, 127, 133, 127, 101, 106, 133, 106, 97, 97, 101, 101, 81, 93, 93, 93, 89, 89, 85, 89, 72, 75, 78, 81, };
static const short divrem_cutoff_tab[] = {166, 139, 139, 139, 69, 75, 89, 139, 127, 111, 111, 127, 116, 111, 106, 101, 97, 93, 106, 106, 85, 81, 78, 101, 54, 75, 85, 81, 63, 60, 58, 75, 52, 58, 58, 66, 52, 58, 52, 56, 54, 54, 46, 54, 52, 50, 46, 50, 39, 46, 48, 46, 46, 38, 44, 44, 40, 40, 42, 44, };

int
_mpn_mod_poly_inv_series(nn_ptr Q, nn_srcptr B, slong lenB, slong len, gr_ctx_t ctx)
{
    slong tab_i, cutoff, bits;

    lenB = FLINT_MIN(len, lenB);

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

    bits = MPN_MOD_CTX_MODULUS_BITS(ctx);
    tab_i = (bits - FLINT_BITS - 1) / 16;
    cutoff = inv_series_cutoff_tab[tab_i];

    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);
}

int
_mpn_mod_poly_div_series(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, slong len, gr_ctx_t ctx)
{
    slong tab_i, 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 = MPN_MOD_CTX_MODULUS_BITS(ctx);
    tab_i = (bits - FLINT_BITS - 1) / 16;
    cutoff = div_series_cutoff_tab[tab_i];

    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);
}

int
_mpn_mod_poly_div(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx)
{
    slong tab_i, cutoff, bits;

    bits = MPN_MOD_CTX_MODULUS_BITS(ctx);
    tab_i = (bits - FLINT_BITS - 1) / 16;
    cutoff = div_series_cutoff_tab[tab_i];

    if (lenB <= cutoff || lenA - lenB + 1 <= cutoff)
        return _gr_poly_div_basecase(Q, A, lenA, B, lenB, ctx);
    else
        return _gr_poly_div_newton(Q, A, lenA, B, lenB, ctx);
}

/* note: we don't define _mpn_mod_poly_divexact because the default algorithm is currently fine */

/* todo: check unbalanced tuning */
int _mpn_mod_poly_divrem(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx)
{
    slong tab_i, cutoff, bits;

    bits = MPN_MOD_CTX_MODULUS_BITS(ctx);
    tab_i = (bits - FLINT_BITS - 1) / 16;
    cutoff = divrem_cutoff_tab[tab_i];

    if (lenA <= cutoff || lenB <= cutoff || lenA - lenB <= cutoff)
        return _mpn_mod_poly_divrem_basecase(Q, R, A, lenA, B, lenB, ctx);
    else
        return _gr_poly_divrem_newton(Q, R, A, lenA, B, lenB, ctx);
}

int _mpn_mod_poly_gcd(nn_ptr G, slong * lenG, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx)
{
    slong cutoff = 240;

    if (lenA < cutoff || lenB < 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, cutoff / 3, cutoff, ctx);
}

int _mpn_mod_poly_xgcd(slong * lenG, nn_ptr G, nn_ptr S, nn_ptr T, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx)
{
    slong cutoff = 240;

    if (lenA < cutoff || lenB < cutoff)
        return _gr_poly_xgcd_euclidean(lenG, G, S, T, A, lenA, B, lenB, ctx);
    else
        return _gr_poly_xgcd_hgcd(lenG, G, S, T, A, lenA, B, lenB, cutoff / 3, cutoff, ctx);
}