From 76900436010fc1cd670ddf676c1f6958340c7d5d Mon Sep 17 00:00:00 2001 From: Ricardo Buring Date: Fri, 8 Sep 2023 13:24:59 +0200 Subject: [PATCH 1/4] Fix memory leak in fmpz_mod_poly_inv_series --- src/fmpz_mod_poly/inv_series.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/fmpz_mod_poly/inv_series.c b/src/fmpz_mod_poly/inv_series.c index 7e0772c342..d85de05651 100644 --- a/src/fmpz_mod_poly/inv_series.c +++ b/src/fmpz_mod_poly/inv_series.c @@ -46,6 +46,7 @@ void fmpz_mod_poly_inv_series(fmpz_mod_poly_t g, const fmpz_mod_poly_t h, slong _fmpz_mod_poly_set_length(t, n); _fmpz_mod_poly_normalise(t); fmpz_mod_poly_swap(g, t, ctx); + fmpz_mod_poly_clear(t, ctx); } else { From fef2cb563d81d863108a0d9479aadd6345384c1b Mon Sep 17 00:00:00 2001 From: Ricardo Buring Date: Fri, 8 Sep 2023 16:59:41 +0200 Subject: [PATCH 2/4] Fix some memory leaks in tests --- src/acb_poly/test/t-inv_series.c | 1 + src/acb_poly/test/t-zeta_em_tail_bsplit.c | 1 + src/aprcl/test/t-unity_zp_reduce_cyclotomic.c | 1 + src/bernoulli/test/t-rev.c | 2 ++ src/fmpz/test/t-xgcd_canonical_bezout.c | 1 + src/fq_nmod/test/t-mul_si.c | 2 ++ src/fq_templates/test/t-sqrt.c | 1 + src/gr_mat/test/t-inv.c | 1 + 8 files changed, 10 insertions(+) diff --git a/src/acb_poly/test/t-inv_series.c b/src/acb_poly/test/t-inv_series.c index 5370bc51b0..2ad3458d2a 100644 --- a/src/acb_poly/test/t-inv_series.c +++ b/src/acb_poly/test/t-inv_series.c @@ -117,6 +117,7 @@ int main(void) acb_poly_clear(a); acb_poly_clear(b); acb_poly_clear(ab); + fmpq_poly_clear(id); } flint_randclear(state); diff --git a/src/acb_poly/test/t-zeta_em_tail_bsplit.c b/src/acb_poly/test/t-zeta_em_tail_bsplit.c index 1d026ce302..76fbad1ba3 100644 --- a/src/acb_poly/test/t-zeta_em_tail_bsplit.c +++ b/src/acb_poly/test/t-zeta_em_tail_bsplit.c @@ -63,6 +63,7 @@ int main(void) acb_clear(Na); acb_clear(s); + _acb_vec_clear(Nasx, len); _acb_vec_clear(z1, len); _acb_vec_clear(z2, len); } diff --git a/src/aprcl/test/t-unity_zp_reduce_cyclotomic.c b/src/aprcl/test/t-unity_zp_reduce_cyclotomic.c index bd76101daa..3f1d063298 100644 --- a/src/aprcl/test/t-unity_zp_reduce_cyclotomic.c +++ b/src/aprcl/test/t-unity_zp_reduce_cyclotomic.c @@ -83,6 +83,7 @@ int main(void) unity_zp_clear(g); } + fmpz_mod_ctx_clear(ctx); FLINT_TEST_CLEANUP(state); flint_printf("PASS\n"); diff --git a/src/bernoulli/test/t-rev.c b/src/bernoulli/test/t-rev.c index 9526cfabfd..b12bc837c5 100644 --- a/src/bernoulli/test/t-rev.c +++ b/src/bernoulli/test/t-rev.c @@ -89,6 +89,8 @@ int main(void) fmpz_clear(denom); } + nmod_poly_clear(A); + flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); diff --git a/src/fmpz/test/t-xgcd_canonical_bezout.c b/src/fmpz/test/t-xgcd_canonical_bezout.c index e28767b92c..176fff6bb9 100644 --- a/src/fmpz/test/t-xgcd_canonical_bezout.c +++ b/src/fmpz/test/t-xgcd_canonical_bezout.c @@ -544,6 +544,7 @@ main(void) fmpz_clear(G); } + fmpz_clear(maxval); fmpz_clear(nd); fmpz_clear(na); fmpz_clear(nb); diff --git a/src/fq_nmod/test/t-mul_si.c b/src/fq_nmod/test/t-mul_si.c index 641a05aa89..25673e9e5e 100644 --- a/src/fq_nmod/test/t-mul_si.c +++ b/src/fq_nmod/test/t-mul_si.c @@ -50,6 +50,8 @@ main(void) } } + fq_nmod_clear(rop, ctx); + fq_nmod_ctx_clear(ctx); fmpz_clear(p); fmpz_clear(f); diff --git a/src/fq_templates/test/t-sqrt.c b/src/fq_templates/test/t-sqrt.c index 0142a0d63f..c9b6401910 100644 --- a/src/fq_templates/test/t-sqrt.c +++ b/src/fq_templates/test/t-sqrt.c @@ -113,6 +113,7 @@ main(void) TEMPLATE(T, clear)(b, ctx); TEMPLATE(T, clear)(c, ctx); TEMPLATE(T, clear)(d, ctx); + TEMPLATE(T, clear)(x, ctx); } TEMPLATE(T, ctx_clear)(ctx); diff --git a/src/gr_mat/test/t-inv.c b/src/gr_mat/test/t-inv.c index ea9a9deabc..07281b7e82 100644 --- a/src/gr_mat/test/t-inv.c +++ b/src/gr_mat/test/t-inv.c @@ -93,6 +93,7 @@ int main(void) gr_mat_clear(A, ctx); gr_mat_clear(B, ctx); + gr_mat_clear(AB, ctx); gr_ctx_clear(ctx); } From 12b0656a3374e79bd4367fff2e384f2043ca4337 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 6 Sep 2023 14:11:38 +0200 Subject: [PATCH 3/4] more divexact algorithms --- doc/source/gr_poly.rst | 8 +- src/gr_poly.h | 12 ++ src/gr_poly/divexact_bidirectional.c | 166 +++++++++++++++++++++++++++ 3 files changed, 185 insertions(+), 1 deletion(-) create mode 100644 src/gr_poly/divexact_bidirectional.c diff --git a/doc/source/gr_poly.rst b/doc/source/gr_poly.rst index 35632ce0a1..1a06930b96 100644 --- a/doc/source/gr_poly.rst +++ b/doc/source/gr_poly.rst @@ -295,10 +295,16 @@ These functions compute a quotient `A / B` which is known to be exact (without remainder) in `R[x]` (or in `R[[x]] / x^n` in the case of series division). Given a nonexact division, they may return gibberish instead of returning an error flag. +For checked exact division, use the ``div`` functions instead. + There are no ``preinv1`` versions because those are identical to the ``div`` counterparts. -.. function:: int _gr_poly_divexact_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx) +.. function:: int _gr_poly_divexact_basecase_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx) + int gr_poly_divexact_basecase_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) + int _gr_poly_divexact_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx) + int gr_poly_divexact_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) + int _gr_poly_divexact_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx) int _gr_poly_divexact_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx) int gr_poly_divexact_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) diff --git a/src/gr_poly.h b/src/gr_poly.h index b0fb47c2ec..646a3c8eba 100644 --- a/src/gr_poly.h +++ b/src/gr_poly.h @@ -227,6 +227,18 @@ GR_POLY_INLINE WARN_UNUSED_RESULT int _gr_poly_div_series(gr_ptr res, gr_srcptr WARN_UNUSED_RESULT int _gr_poly_div_series_generic(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx); WARN_UNUSED_RESULT int gr_poly_div_series(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_poly_divexact_basecase_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_poly_divexact_basecase_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_poly_divexact_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_poly_divexact_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_poly_divexact_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_poly_divexact_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_poly_divexact_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); + +WARN_UNUSED_RESULT int _gr_poly_divexact_series_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx); +WARN_UNUSED_RESULT int _gr_poly_divexact_series_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_poly_divexact_series_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx); + WARN_UNUSED_RESULT int _gr_poly_sqrt_series_newton(gr_ptr res, gr_srcptr f, slong flen, slong len, slong cutoff, gr_ctx_t ctx); WARN_UNUSED_RESULT int gr_poly_sqrt_series_newton(gr_poly_t res, const gr_poly_t f, slong len, slong cutoff, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_poly_sqrt_series_basecase(gr_ptr res, gr_srcptr f, slong flen, slong len, gr_ctx_t ctx); diff --git a/src/gr_poly/divexact_bidirectional.c b/src/gr_poly/divexact_bidirectional.c new file mode 100644 index 0000000000..55722b4ccd --- /dev/null +++ b/src/gr_poly/divexact_bidirectional.c @@ -0,0 +1,166 @@ +/* + 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 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "gr_vec.h" +#include "gr_poly.h" + +static int +__gr_poly_divexact_bidirectional(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, int basecase, gr_ctx_t ctx) +{ + slong lenQ, len_lo, len_hi; + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + + /* See if we can divide from the bottom. */ + while (lenB > 1) + { + truth_t is_zero; + + is_zero = gr_is_zero(B, ctx); + + if (is_zero == T_FALSE) + break; + + /* We cannot tell if the low coefficient of B is zero, + so series division will not work. However, dividing + from the top might still be possible. */ + if (is_zero == T_UNKNOWN) + { + if (basecase) + return _gr_poly_divexact_basecase(Q, A, lenA, B, lenB, ctx); + else + return _gr_poly_div(Q, A, lenA, B, lenB, ctx); + } + + /* Discard trailing zero coefficient. B is assumed to divide A, + there is a corresponding zero coefficient in A, + which we don't need to check since the contract for divexact + allows returning nonsense in case of an inexact division. */ + B = GR_ENTRY(B, 1, sz); + A = GR_ENTRY(A, 1, sz); + lenB--; + lenA--; + } + + if (lenB == 1) + { + return _gr_vec_divexact_scalar(Q, A, lenA, B, ctx); + } + + lenQ = lenA - lenB + 1; + len_hi = lenQ / 2; + len_lo = lenQ - len_hi; + + if (basecase) + { + status |= _gr_poly_div_series(Q, A, lenA, B, lenB, len_lo, ctx); + status |= _gr_poly_div(GR_ENTRY(Q, len_lo, sz), GR_ENTRY(A, len_lo, sz), lenA - len_lo, B, lenB, ctx); + } + else + { + status |= _gr_poly_divexact_series_basecase(Q, A, lenA, B, lenB, len_lo, ctx); + status |= _gr_poly_divexact_basecase(GR_ENTRY(Q, len_lo, sz), GR_ENTRY(A, len_lo, sz), lenA - len_lo, B, lenB, ctx); + } + + return status; +} + +int +_gr_poly_divexact_bidirectional(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) +{ + return __gr_poly_divexact_bidirectional(Q, A, lenA, B, lenB, 0, ctx); +} + +int +_gr_poly_divexact_basecase_bidirectional(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) +{ + return __gr_poly_divexact_bidirectional(Q, A, lenA, B, lenB, 1, ctx); +} + +int +gr_poly_divexact_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) +{ + slong Alen, Blen, Qlen; + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + + Alen = A->length; + Blen = B->length; + + if (Blen == 0) + return GR_DOMAIN; + + if (gr_is_zero(GR_ENTRY(B->coeffs, Blen - 1, sz), ctx) != T_FALSE) + return GR_UNABLE; + + if (Alen < Blen) + return gr_poly_zero(Q, ctx); + + Qlen = Alen - Blen + 1; + + if (Q == A || Q == B) + { + gr_poly_t t; + gr_poly_init2(t, Qlen, ctx); + status = _gr_poly_divexact_bidirectional(t->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx); + gr_poly_swap(Q, t, ctx); + gr_poly_clear(t, ctx); + } + else + { + gr_poly_fit_length(Q, Qlen, ctx); + status = _gr_poly_divexact_bidirectional(Q->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx); + } + + _gr_poly_set_length(Q, Qlen, ctx); + _gr_poly_normalise(Q, ctx); + return status; +} + +int +gr_poly_divexact_basecase_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) +{ + slong Alen, Blen, Qlen; + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + + Alen = A->length; + Blen = B->length; + + if (Blen == 0) + return GR_DOMAIN; + + if (gr_is_zero(GR_ENTRY(B->coeffs, Blen - 1, sz), ctx) != T_FALSE) + return GR_UNABLE; + + if (Alen < Blen) + return gr_poly_zero(Q, ctx); + + Qlen = Alen - Blen + 1; + + if (Q == A || Q == B) + { + gr_poly_t t; + gr_poly_init2(t, Qlen, ctx); + status = _gr_poly_divexact_bidirectional(t->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx); + gr_poly_swap(Q, t, ctx); + gr_poly_clear(t, ctx); + } + else + { + gr_poly_fit_length(Q, Qlen, ctx); + status = _gr_poly_divexact_bidirectional(Q->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx); + } + + _gr_poly_set_length(Q, Qlen, ctx); + _gr_poly_normalise(Q, ctx); + return status; +} From 58d3aacc78558a35391d7af76822b128669c8e88 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 6 Sep 2023 21:19:45 +0200 Subject: [PATCH 4/4] more gr_poly division work --- doc/source/gr_poly.rst | 35 +++++--- src/gr.h | 2 + src/gr_generic/generic.c | 2 + src/gr_poly.h | 9 +- src/gr_poly/div.c | 2 +- src/gr_poly/divexact.c | 62 +++++++++++++ src/gr_poly/divexact_bidirectional.c | 8 +- src/gr_poly/test/t-divexact.c | 126 +++++++++++++++++++++++++++ 8 files changed, 230 insertions(+), 16 deletions(-) create mode 100644 src/gr_poly/divexact.c create mode 100644 src/gr_poly/test/t-divexact.c diff --git a/doc/source/gr_poly.rst b/doc/source/gr_poly.rst index 1a06930b96..cc340a3fc5 100644 --- a/doc/source/gr_poly.rst +++ b/doc/source/gr_poly.rst @@ -277,8 +277,8 @@ Division uses the Karp-Markstein algorithm. .. function:: int _gr_poly_div_series_newton(gr_ptr res, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, slong cutoff, gr_ctx_t ctx) int gr_poly_div_series_newton(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx) - int _gr_poly_div_series_divconquer(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx); - int gr_poly_div_series_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx); + int _gr_poly_div_series_divconquer(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, slong cutoff, gr_ctx_t ctx) + int gr_poly_div_series_divconquer(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, slong cutoff, gr_ctx_t ctx) int _gr_poly_div_series_invmul(gr_ptr res, gr_srcptr B, slong Blen, gr_srcptr A, slong Alen, slong len, gr_ctx_t ctx) int gr_poly_div_series_invmul(gr_poly_t res, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx) int _gr_poly_div_series_basecase_preinv1(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_srcptr Binv, slong len, gr_ctx_t ctx) @@ -291,14 +291,29 @@ Division uses the Karp-Markstein algorithm. Exact division -------------------------------------------------------------------------------- -These functions compute a quotient `A / B` which is known to be exact +These functions compute a quotient `Q = A / B` which is known to be exact (without remainder) in `R[x]` (or in `R[[x]] / x^n` in the case of series -division). Given a nonexact division, they may return gibberish instead of -returning an error flag. -For checked exact division, use the ``div`` functions instead. +division). Given a nonexact division, they are allowed to set `Q` to +an arbitrary polynomial and return ``GR_SUCCESS`` instead of returning an +error flag. -There are no ``preinv1`` versions because those are identical to the -``div`` counterparts. +`R` is assumed to be an integral domain (this is not checked). + +For exact division, we have the choice of starting the division +from the most significant terms (classical division) or the least significant +(power series division). Which direction is more efficient depends +in part on whether the leading or trailing coefficient of `B` is cheaper +to use for divisions. In a generic setting, this is hard to predict. + +The *bidirectional* algorithms combine two half-divisions from both ends. +This halves the number of operations in the basecase regime, though an +extra coefficient inversion may be needed. + +The ``noinv`` versions perform repeated ``divexact`` operations in the +scalar domain without attempting to invert the leading (or trailing) coefficient, +while other versions check invertibility first. +There are no ``divexact_preinv1`` versions because those are identical to the +``div_preinv1`` counterparts. .. function:: int _gr_poly_divexact_basecase_bidirectional(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx) int gr_poly_divexact_basecase_bidirectional(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) @@ -500,8 +515,8 @@ GCD .. function:: int _gr_poly_xgcd_euclidean(slong * lenG, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_srcptr invB, gr_ctx_t ctx) int gr_poly_xgcd_euclidean(gr_poly_t G, gr_poly_t S, gr_poly_t T, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) -.. function:: int _gr_poly_xgcd_hgcd(slong * Glen, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong hgcd_cutoff, slong cutoff, gr_ctx_t ctx); - int gr_poly_xgcd_hgcd(gr_poly_t G, gr_poly_t S, gr_poly_t T, const gr_poly_t A, const gr_poly_t B, slong hgcd_cutoff, slong cutoff, gr_ctx_t ctx); +.. function:: int _gr_poly_xgcd_hgcd(slong * Glen, gr_ptr G, gr_ptr S, gr_ptr T, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, slong hgcd_cutoff, slong cutoff, gr_ctx_t ctx) + int gr_poly_xgcd_hgcd(gr_poly_t G, gr_poly_t S, gr_poly_t T, const gr_poly_t A, const gr_poly_t B, slong hgcd_cutoff, slong cutoff, gr_ctx_t ctx) Resultant ------------------------------------------------------------------------------- diff --git a/src/gr.h b/src/gr.h index 4087e178b5..c3a96b19e0 100644 --- a/src/gr.h +++ b/src/gr.h @@ -663,7 +663,9 @@ typedef enum /* Polynomial methods (todo: rename -> GR_POLY) */ GR_METHOD_POLY_MULLOW, + GR_METHOD_POLY_DIV, GR_METHOD_POLY_DIVREM, + GR_METHOD_POLY_DIVEXACT, GR_METHOD_POLY_TAYLOR_SHIFT, GR_METHOD_POLY_INV_SERIES, GR_METHOD_POLY_INV_SERIES_BASECASE, diff --git a/src/gr_generic/generic.c b/src/gr_generic/generic.c index 079b06a6a7..c1dfe92680 100644 --- a/src/gr_generic/generic.c +++ b/src/gr_generic/generic.c @@ -2702,7 +2702,9 @@ const gr_method_tab_input _gr_generic_methods[] = {GR_METHOD_VEC_SET_POWERS, (gr_funcptr) gr_generic_vec_set_powers}, {GR_METHOD_POLY_MULLOW, (gr_funcptr) _gr_poly_mullow_generic}, + {GR_METHOD_POLY_DIV, (gr_funcptr) _gr_poly_div_generic}, {GR_METHOD_POLY_DIVREM, (gr_funcptr) _gr_poly_divrem_generic}, + {GR_METHOD_POLY_DIVEXACT, (gr_funcptr) _gr_poly_divexact_generic}, {GR_METHOD_POLY_TAYLOR_SHIFT, (gr_funcptr) _gr_poly_taylor_shift_generic}, {GR_METHOD_POLY_INV_SERIES, (gr_funcptr) _gr_poly_inv_series_generic}, {GR_METHOD_POLY_INV_SERIES_BASECASE,(gr_funcptr) _gr_poly_inv_series_basecase_generic}, diff --git a/src/gr_poly.h b/src/gr_poly.h index 646a3c8eba..7a5994f942 100644 --- a/src/gr_poly.h +++ b/src/gr_poly.h @@ -189,7 +189,9 @@ WARN_UNUSED_RESULT int gr_poly_div_basecase(gr_poly_t Q, const gr_poly_t A, cons WARN_UNUSED_RESULT int _gr_poly_div_newton(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx); WARN_UNUSED_RESULT int gr_poly_div_newton(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); -WARN_UNUSED_RESULT int _gr_poly_div(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx); + +GR_POLY_INLINE WARN_UNUSED_RESULT int _gr_poly_div(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) { return GR_POLY_BINARY_OP(ctx, POLY_DIV)(Q, A, lenA, B, lenB, ctx); } +WARN_UNUSED_RESULT int _gr_poly_div_generic(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx); WARN_UNUSED_RESULT int gr_poly_div(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_poly_rem(gr_ptr R, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx); @@ -235,6 +237,11 @@ WARN_UNUSED_RESULT int _gr_poly_divexact_basecase_noinv(gr_ptr Q, gr_srcptr A, s WARN_UNUSED_RESULT int _gr_poly_divexact_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, gr_ctx_t ctx); WARN_UNUSED_RESULT int gr_poly_divexact_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); +GR_POLY_INLINE WARN_UNUSED_RESULT int _gr_poly_divexact(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) { return GR_POLY_BINARY_OP(ctx, POLY_DIVEXACT)(Q, A, lenA, B, lenB, ctx); } +WARN_UNUSED_RESULT int _gr_poly_divexact_generic(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx); +WARN_UNUSED_RESULT int gr_poly_divexact(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx); + + WARN_UNUSED_RESULT int _gr_poly_divexact_series_basecase_noinv(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx); WARN_UNUSED_RESULT int _gr_poly_divexact_series_basecase(gr_ptr Q, gr_srcptr A, slong Alen, gr_srcptr B, slong Blen, slong len, gr_ctx_t ctx); WARN_UNUSED_RESULT int gr_poly_divexact_series_basecase(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, slong len, gr_ctx_t ctx); diff --git a/src/gr_poly/div.c b/src/gr_poly/div.c index 51f230664c..8088232d32 100644 --- a/src/gr_poly/div.c +++ b/src/gr_poly/div.c @@ -17,7 +17,7 @@ #define DIVCONQUER_CUTOFF 10 int -_gr_poly_div(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) +_gr_poly_div_generic(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) { int status; diff --git a/src/gr_poly/divexact.c b/src/gr_poly/divexact.c new file mode 100644 index 0000000000..ad972468cb --- /dev/null +++ b/src/gr_poly/divexact.c @@ -0,0 +1,62 @@ +/* + 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 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "gr_vec.h" +#include "gr_poly.h" + +int +_gr_poly_divexact_generic(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, slong lenB, gr_ctx_t ctx) +{ + if (lenB <= 2) + return _gr_poly_divexact_basecase(Q, A, lenA, B, lenB, ctx); + else + return _gr_poly_divexact_bidirectional(Q, A, lenA, B, lenB, ctx); +} + +int +gr_poly_divexact(gr_poly_t Q, const gr_poly_t A, const gr_poly_t B, gr_ctx_t ctx) +{ + slong Alen, Blen, Qlen; + int status = GR_SUCCESS; + slong sz = ctx->sizeof_elem; + + Alen = A->length; + Blen = B->length; + + if (Blen == 0) + return GR_DOMAIN; + + if (gr_is_zero(GR_ENTRY(B->coeffs, Blen - 1, sz), ctx) != T_FALSE) + return GR_UNABLE; + + if (Alen < Blen) + return gr_poly_zero(Q, ctx); + + Qlen = Alen - Blen + 1; + + if (Q == A || Q == B) + { + gr_poly_t t; + gr_poly_init2(t, Qlen, ctx); + status = _gr_poly_divexact(t->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx); + gr_poly_swap(Q, t, ctx); + gr_poly_clear(t, ctx); + } + else + { + gr_poly_fit_length(Q, Qlen, ctx); + status = _gr_poly_divexact(Q->coeffs, A->coeffs, A->length, B->coeffs, B->length, ctx); + } + + _gr_poly_set_length(Q, Qlen, ctx); + _gr_poly_normalise(Q, ctx); + return status; +} diff --git a/src/gr_poly/divexact_bidirectional.c b/src/gr_poly/divexact_bidirectional.c index 55722b4ccd..672c439a04 100644 --- a/src/gr_poly/divexact_bidirectional.c +++ b/src/gr_poly/divexact_bidirectional.c @@ -61,13 +61,13 @@ __gr_poly_divexact_bidirectional(gr_ptr Q, gr_srcptr A, slong lenA, gr_srcptr B, if (basecase) { - status |= _gr_poly_div_series(Q, A, lenA, B, lenB, len_lo, ctx); - status |= _gr_poly_div(GR_ENTRY(Q, len_lo, sz), GR_ENTRY(A, len_lo, sz), lenA - len_lo, B, lenB, ctx); + status |= _gr_poly_divexact_series_basecase(Q, A, lenA, B, lenB, len_lo, ctx); + status |= _gr_poly_divexact_basecase(GR_ENTRY(Q, len_lo, sz), GR_ENTRY(A, len_lo, sz), lenA - len_lo, B, lenB, ctx); } else { - status |= _gr_poly_divexact_series_basecase(Q, A, lenA, B, lenB, len_lo, ctx); - status |= _gr_poly_divexact_basecase(GR_ENTRY(Q, len_lo, sz), GR_ENTRY(A, len_lo, sz), lenA - len_lo, B, lenB, ctx); + status |= _gr_poly_div_series(Q, A, lenA, B, lenB, len_lo, ctx); + status |= _gr_poly_div(GR_ENTRY(Q, len_lo, sz), GR_ENTRY(A, len_lo, sz), lenA - len_lo, B, lenB, ctx); } return status; diff --git a/src/gr_poly/test/t-divexact.c b/src/gr_poly/test/t-divexact.c new file mode 100644 index 0000000000..78e0b66945 --- /dev/null +++ b/src/gr_poly/test/t-divexact.c @@ -0,0 +1,126 @@ +/* + 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 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "ulong_extras.h" +#include "gr_poly.h" + +int main(void) +{ + slong iter; + flint_rand_t state; + + flint_printf("divexact...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + int status; + gr_ctx_t ctx; + gr_poly_t A, B, Q, Q2; + + gr_ctx_init_random(ctx, state); + + gr_poly_init(A, ctx); + gr_poly_init(B, ctx); + gr_poly_init(Q, ctx); + gr_poly_init(Q2, ctx); + + status = GR_SUCCESS; + + status |= gr_poly_randtest(B, state, 1 + n_randint(state, 6), ctx); + + if (gr_poly_is_zero(B, ctx) == T_FALSE) + { + status |= gr_poly_randtest(Q, state, 1 + n_randint(state, 6), ctx); + status |= gr_poly_randtest(Q2, state, 1 + n_randint(state, 6), ctx); + status |= gr_poly_mul(A, Q2, B, ctx); + + switch (n_randint(state, 12)) + { + case 0: + status |= gr_poly_set(Q, A, ctx); + status |= gr_poly_divexact(Q, Q, B, ctx); + break; + case 1: + status |= gr_poly_set(Q, B, ctx); + status |= gr_poly_divexact(Q, A, B, ctx); + break; + case 2: + status |= gr_poly_divexact(Q, A, B, ctx); + break; + + case 3: + status |= gr_poly_set(Q, A, ctx); + status |= gr_poly_divexact_basecase(Q, Q, B, ctx); + break; + case 4: + status |= gr_poly_set(Q, B, ctx); + status |= gr_poly_divexact_basecase(Q, A, B, ctx); + break; + case 5: + status |= gr_poly_divexact_basecase(Q, A, B, ctx); + break; + + case 6: + status |= gr_poly_set(Q, A, ctx); + status |= gr_poly_divexact_bidirectional(Q, Q, B, ctx); + break; + case 7: + status |= gr_poly_set(Q, B, ctx); + status |= gr_poly_divexact_bidirectional(Q, A, B, ctx); + break; + case 8: + status |= gr_poly_divexact_bidirectional(Q, A, B, ctx); + break; + + case 9: + status |= gr_poly_set(Q, A, ctx); + status |= gr_poly_divexact_basecase_bidirectional(Q, Q, B, ctx); + break; + case 10: + status |= gr_poly_set(Q, B, ctx); + status |= gr_poly_divexact_basecase_bidirectional(Q, A, B, ctx); + break; + default: + status |= gr_poly_divexact_basecase_bidirectional(Q, A, B, ctx); + break; + } + + if (gr_ctx_is_integral_domain(ctx) == T_TRUE && + ((status == GR_SUCCESS && gr_poly_equal(Q, Q2, ctx) == T_FALSE) + || status == GR_DOMAIN + || (ctx->which_ring == GR_CTX_FMPZ && status != GR_SUCCESS))) + { + flint_printf("FAIL\n\n"); + gr_ctx_println(ctx); + flint_printf("A = "); gr_poly_print(A, ctx); flint_printf("\n"); + flint_printf("B = "); gr_poly_print(B, ctx); flint_printf("\n"); + flint_printf("Q = "); gr_poly_print(Q, ctx); flint_printf("\n"); + flint_printf("Q2 = "); gr_poly_print(Q2, ctx); flint_printf("\n"); + flint_abort(); + } + } + + gr_poly_clear(A, ctx); + gr_poly_clear(B, ctx); + gr_poly_clear(Q, ctx); + gr_poly_clear(Q2, ctx); + + gr_ctx_clear(ctx); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return 0; +}