Skip to content

Commit

Permalink
more gr_poly division work
Browse files Browse the repository at this point in the history
  • Loading branch information
fredrik-johansson committed Sep 6, 2023
1 parent c4f6ae2 commit ae81eb6
Show file tree
Hide file tree
Showing 8 changed files with 230 additions and 16 deletions.
35 changes: 25 additions & 10 deletions doc/source/gr_poly.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
-------------------------------------------------------------------------------
Expand Down
2 changes: 2 additions & 0 deletions src/gr.h
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
2 changes: 2 additions & 0 deletions src/gr_generic/generic.c
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand Down
9 changes: 8 additions & 1 deletion src/gr_poly.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion src/gr_poly/div.c
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
62 changes: 62 additions & 0 deletions src/gr_poly/divexact.c
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#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;
}
8 changes: 4 additions & 4 deletions src/gr_poly/divexact_bidirectional.c
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
126 changes: 126 additions & 0 deletions src/gr_poly/test/t-divexact.c
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

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

0 comments on commit ae81eb6

Please sign in to comment.