diff --git a/arith/number_of_partitions_vec.c b/arith/number_of_partitions_vec.c index 7e4e448cbe..0f4ee0ccd8 100644 --- a/arith/number_of_partitions_vec.c +++ b/arith/number_of_partitions_vec.c @@ -51,7 +51,7 @@ arith_number_of_partitions_vec(fmpz * res, slong len) if (n + k < len) tmp[n + k] = WORD(-1); if (n + 3*k + 1 < len) tmp[n + 3*k + 1] = WORD(1); - _fmpz_poly_inv_series(res, tmp, len); + _fmpz_poly_inv_series(res, tmp, len, len); _fmpz_vec_clear(tmp, len); } diff --git a/flint.h b/flint.h index cf64cccf03..1afebe29e4 100644 --- a/flint.h +++ b/flint.h @@ -319,6 +319,30 @@ void mpn_tdiv_q(mp_ptr qp, __tmp_root = __tmp_root->next; \ } +/* Newton iteration macros */ +#define FLINT_NEWTON_INIT(from, to) \ + { \ + slong __steps[FLINT_BITS], __i, __from, __to; \ + __steps[__i = 0] = __to = (to); \ + __from = (from); \ + while (__to > __from) \ + __steps[++__i] = (__to = (__to + 1) / 2); \ + +#define FLINT_NEWTON_BASECASE(bc_to) { slong bc_to = __to; + +#define FLINT_NEWTON_END_BASECASE } + +#define FLINT_NEWTON_LOOP(step_from, step_to) \ + { \ + for (__i--; __i >= 0; __i--) \ + { \ + slong step_from = __steps[__i+1]; \ + slong step_to = __steps[__i]; \ + +#define FLINT_NEWTON_END_LOOP }} + +#define FLINT_NEWTON_END } + int parse_fmt(int * floating, const char * fmt); size_t flint_printf(const char * str, ...); /* flint version of printf */ diff --git a/fmpz_poly.h b/fmpz_poly.h index b896249650..11a2975c18 100644 --- a/fmpz_poly.h +++ b/fmpz_poly.h @@ -47,7 +47,7 @@ extern "C" { #endif -#define FMPZ_POLY_INV_NEWTON_CUTOFF 32 +#define FMPZ_POLY_INV_NEWTON_CUTOFF 32 /* Type definitions *********************************************************/ @@ -717,21 +717,17 @@ _fmpz_poly_div_root(fmpz * Q, const fmpz * A, slong len, const fmpz_t c); /* Power series division ***************************************************/ -void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n); +void _fmpz_poly_inv_series_basecase(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n); + +void fmpz_poly_inv_series_basecase(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n); + +void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n); void fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n); -static __inline__ void -_fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong n) -{ - _fmpz_poly_inv_series_newton(Qinv, Q, n); -} +void _fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n); -static __inline__ void -fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) -{ - fmpz_poly_inv_series_newton(Qinv, Q, n); -} +void fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n); void _fmpz_poly_div_series(fmpz * Q, const fmpz * A, const fmpz * B, slong n); diff --git a/fmpz_poly/div_series.c b/fmpz_poly/div_series.c index a84c8b19f0..53f9248486 100644 --- a/fmpz_poly/div_series.c +++ b/fmpz_poly/div_series.c @@ -41,7 +41,7 @@ _fmpz_poly_div_series(fmpz * Q, const fmpz * A, const fmpz * B, slong n) { fmpz * Binv = _fmpz_vec_init(n); - _fmpz_poly_inv_series(Binv, B, n); + _fmpz_poly_inv_series(Binv, B, n, n); _fmpz_poly_mullow(Q, A, n, Binv, n, n); _fmpz_vec_clear(Binv, n); diff --git a/fmpz_poly/doc/fmpz_poly.txt b/fmpz_poly/doc/fmpz_poly.txt index 40510018a8..4cac74dc38 100644 --- a/fmpz_poly/doc/fmpz_poly.txt +++ b/fmpz_poly/doc/fmpz_poly.txt @@ -1694,31 +1694,45 @@ int fmpz_poly_divides(fmpz_poly_t Q, ******************************************************************************* -void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n) +void _fmpz_poly_inv_series_basecase(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n) - Computes the first $n$ terms of the inverse power series of $Q$ using - Newton iteration. + Computes the first $n$ terms of the inverse power series of + \code{(Q, lenQ)} using a recurrence. - Assumes that $n \geq 1$, that $Q$ has length at least $n$ and constant - term~$\pm 1$. Does not support aliasing. + Assumes that $n \geq 1$ and that $Q$ has constant term~$\pm 1$. + Does not support aliasing. -id fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) +void fmpz_poly_inv_series_basecase(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) Computes the first $n$ terms of the inverse power series of $Q$ - using Newton iteration, assuming that $Q$ has constant term~$\pm 1$ + using a recurrence, assuming that $Q$ has constant term~$\pm 1$ and $n \geq 1$. +void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n) + + Computes the first $n$ terms of the inverse power series of + \code{(Q, lenQ)} using Newton iteration. + + Assumes that $n \geq 1$ and that $Q$ has constant term~$\pm 1$. + Does not support aliasing. + +void fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong Qlen, slong n) + + Computes the first $n$ terms of the inverse power series of $Q$ using + Newton iteration, assuming $Q$ has constant term~$\pm 1$ and $n \geq 1$. + void _fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong n) - Computes the first $n$ terms of the inverse power series of $Q$. + Computes the first $n$ terms of the inverse power series of + \code{(Q, lenQ)}. - Assumes that $n \geq 1$, that $Q$ has length at least $n$ and constant - term~$1$. Does not support aliasing. + Assumes that $n \geq 1$ and that $Q$ has constant term~$\pm 1$. + Does not support aliasing. void fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) Computes the first $n$ terms of the inverse power series of $Q$, - assuming $Q$ has constant term~$1$ and $n \geq 1$. + assuming $Q$ has constant term~$\pm 1$ and $n \geq 1$. void _fmpz_poly_div_series(fmpz * Q, const fmpz * A, const fmpz * B) diff --git a/fmpz_poly/eta_qexp.c b/fmpz_poly/eta_qexp.c index b51ab18453..0f63f835ab 100644 --- a/fmpz_poly/eta_qexp.c +++ b/fmpz_poly/eta_qexp.c @@ -176,7 +176,7 @@ _fmpz_poly_eta_qexp(fmpz * f, slong e, slong len) { fmpz * t = _fmpz_vec_init(len); _fmpz_poly_eta_qexp(t, -e, len); - _fmpz_poly_inv_series(f, t, len); + _fmpz_poly_inv_series(f, t, len, len); _fmpz_vec_clear(t, len); return; } diff --git a/fmpz_poly/inv_series.c b/fmpz_poly/inv_series.c new file mode 100644 index 0000000000..f3c445b3a5 --- /dev/null +++ b/fmpz_poly/inv_series.c @@ -0,0 +1,67 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + FLINT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with FLINT; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "fmpz_poly.h" + +void +_fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n) +{ + if (Qlen <= 8 || n <= 24) + _fmpz_poly_inv_series_basecase(Qinv, Q, Qlen, n); + else + _fmpz_poly_inv_series_newton(Qinv, Q, Qlen, n); +} + +void +fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) +{ + slong Qlen = Q->length; + + Qlen = FLINT_MIN(Qlen, n); + + if (Qlen == 0) + { + flint_printf("Exception (fmpz_poly_inv_series). Division by zero.\n"); + abort(); + } + + if (Qinv != Q) + { + fmpz_poly_fit_length(Qinv, n); + _fmpz_poly_inv_series(Qinv->coeffs, Q->coeffs, Qlen, n); + } + else + { + fmpz_poly_t t; + fmpz_poly_init2(t, n); + _fmpz_poly_inv_series(t->coeffs, Q->coeffs, Qlen, n); + fmpz_poly_swap(Qinv, t); + fmpz_poly_clear(t); + } + + _fmpz_poly_set_length(Qinv, n); + _fmpz_poly_normalise(Qinv); +} + diff --git a/fmpz_poly/inv_series_basecase.c b/fmpz_poly/inv_series_basecase.c new file mode 100644 index 0000000000..0d2c42a0f1 --- /dev/null +++ b/fmpz_poly/inv_series_basecase.c @@ -0,0 +1,86 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + FLINT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with FLINT; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "fmpz_poly.h" + +void +_fmpz_poly_inv_series_basecase(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n) +{ + Qlen = FLINT_MIN(Qlen, n); + + fmpz_set(Qinv, Q); + + if (Qlen == 1) + { + _fmpz_vec_zero(Qinv + 1, n - 1); + } + else + { + slong i, j; + + for (i = 1; i < n; i++) + { + fmpz_mul(Qinv + i, Q + 1, Qinv + i - 1); + + for (j = 2; j < FLINT_MIN(i + 1, Qlen); j++) + fmpz_addmul(Qinv + i, Q + j, Qinv + i - j); + + if (fmpz_is_one(Qinv)) + fmpz_neg(Qinv + i, Qinv + i); + } + } +} + +void +fmpz_poly_inv_series_basecase(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) +{ + slong Qlen = Q->length; + + Qlen = FLINT_MIN(Qlen, n); + + if (Qlen == 0) + { + flint_printf("Exception (fmpz_poly_inv_series_basecase). Division by zero.\n"); + abort(); + } + + if (Qinv != Q) + { + fmpz_poly_fit_length(Qinv, n); + _fmpz_poly_inv_series_basecase(Qinv->coeffs, Q->coeffs, Qlen, n); + } + else + { + fmpz_poly_t t; + fmpz_poly_init2(t, n); + _fmpz_poly_inv_series_basecase(t->coeffs, Q->coeffs, Qlen, n); + fmpz_poly_swap(Qinv, t); + fmpz_poly_clear(t); + } + + _fmpz_poly_set_length(Qinv, n); + _fmpz_poly_normalise(Qinv); +} + diff --git a/fmpz_poly/inv_series_newton.c b/fmpz_poly/inv_series_newton.c index 2002062ff4..c639dad947 100644 --- a/fmpz_poly/inv_series_newton.c +++ b/fmpz_poly/inv_series_newton.c @@ -20,102 +20,104 @@ /****************************************************************************** Copyright (C) 2010, 2011 Sebastian Pancratz + Copyright (C) 2014 Fredrik Johansson ******************************************************************************/ -#include -#include -#include "flint.h" -#include "fmpz.h" -#include "fmpz_vec.h" #include "fmpz_poly.h" -void -_fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n) +/* Requires 2*min(Qlen,n) + n - 1 < 3n coefficients of scratch space in W */ +static void +_fmpz_poly_inv_series_basecase_rev(fmpz * Qinv, fmpz * W, + const fmpz * Q, slong Qlen, slong n) { - if (n == 1) /* Q is +-1 */ + slong Wlen; + fmpz *Qrev; + + Qlen = FLINT_MIN(Qlen, n); + Wlen = n + Qlen - 1; + Qrev = W + Wlen; + + _fmpz_poly_reverse(Qrev, Q, Qlen, Qlen); + _fmpz_vec_zero(W, Wlen - 1); + fmpz_one(W + Wlen - 1); + _fmpz_poly_div_basecase(Qinv, W, W, Wlen, Qrev, Qlen); + _fmpz_poly_reverse(Qinv, Qinv, n, n); +} + +#define MULLOW(z, x, xn, y, yn, nn) \ + if ((xn) >= (yn)) \ + _fmpz_poly_mullow(z, x, xn, y, yn, nn); \ + else \ + _fmpz_poly_mullow(z, y, yn, x, xn, nn); \ + +void +_fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n) +{ + Qlen = FLINT_MIN(Qlen, n); + + if (Qlen == 1) { fmpz_set(Qinv, Q); + _fmpz_vec_zero(Qinv + 1, n - 1); } else { - const slong alloc = FLINT_MAX(n, 3 * FMPZ_POLY_INV_NEWTON_CUTOFF); - slong *a, i, m; - fmpz *W; + slong alloc, Qnlen, Wlen, W2len; + fmpz * W; + alloc = FLINT_MAX(n, 3 * FMPZ_POLY_INV_NEWTON_CUTOFF); W = _fmpz_vec_init(alloc); - for (i = 1; (WORD(1) << i) < n; i++) ; - - a = (slong *) flint_malloc(i * sizeof(slong)); - a[i = 0] = n; - while (n >= FMPZ_POLY_INV_NEWTON_CUTOFF) - a[++i] = (n = (n + 1) / 2); - - /* Base case */ - { - fmpz *Qrev = W + 2 * FMPZ_POLY_INV_NEWTON_CUTOFF; - - _fmpz_poly_reverse(Qrev, Q, n, n); - _fmpz_vec_zero(W, 2*n - 2); - fmpz_one(W + (2*n - 2)); - _fmpz_poly_div_basecase(Qinv, W, W, 2*n - 1, Qrev, n); - _fmpz_poly_reverse(Qinv, Qinv, n, n); - } - - for (i--; i >= 0; i--) - { - m = n; - n = a[i]; - - _fmpz_poly_mullow(W, Q, n, Qinv, m, n); - _fmpz_poly_mullow(Qinv + m, Qinv, m, W + m, n - m, n - m); - _fmpz_vec_neg(Qinv + m, Qinv + m, n - m); - } + FLINT_NEWTON_INIT(FMPZ_POLY_INV_NEWTON_CUTOFF, n) + + FLINT_NEWTON_BASECASE(n) + _fmpz_poly_inv_series_basecase_rev(Qinv, W, Q, Qlen, n); + FLINT_NEWTON_END_BASECASE + + FLINT_NEWTON_LOOP(m, n) + Qnlen = FLINT_MIN(Qlen, n); + Wlen = FLINT_MIN(Qnlen + m - 1, n); + W2len = Wlen - m; + MULLOW(W, Q, Qnlen, Qinv, m, Wlen); + MULLOW(Qinv + m, Qinv, m, W + m, W2len, n - m); + _fmpz_vec_neg(Qinv + m, Qinv + m, n - m); + FLINT_NEWTON_END_LOOP + + FLINT_NEWTON_END _fmpz_vec_clear(W, alloc); - flint_free(a); } } -void fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) +void +fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n) { - fmpz *Qcopy; - int Qalloc; + slong Qlen = Q->length; - if (Q->length >= n) - { - Qcopy = Q->coeffs; - Qalloc = 0; - } - else + Qlen = FLINT_MIN(Qlen, n); + + if (Qlen == 0) { - slong i; - Qcopy = (fmpz *) flint_malloc(n * sizeof(fmpz)); - for (i = 0; i < Q->length; i++) - Qcopy[i] = Q->coeffs[i]; - flint_mpn_zero((mp_ptr) Qcopy + i, n - i); - Qalloc = 1; + flint_printf("Exception (fmpz_poly_inv_series_newton). Division by zero.\n"); + abort(); } if (Qinv != Q) { fmpz_poly_fit_length(Qinv, n); - _fmpz_poly_inv_series_newton(Qinv->coeffs, Qcopy, n); + _fmpz_poly_inv_series_newton(Qinv->coeffs, Q->coeffs, Qlen, n); } else { fmpz_poly_t t; fmpz_poly_init2(t, n); - _fmpz_poly_inv_series_newton(t->coeffs, Qcopy, n); + _fmpz_poly_inv_series_newton(t->coeffs, Q->coeffs, Qlen, n); fmpz_poly_swap(Qinv, t); fmpz_poly_clear(t); } - + _fmpz_poly_set_length(Qinv, n); _fmpz_poly_normalise(Qinv); - - if (Qalloc) - flint_free(Qcopy); } diff --git a/fmpz_poly/revert_series_lagrange.c b/fmpz_poly/revert_series_lagrange.c index 886f55acda..e4afeb7cf5 100644 --- a/fmpz_poly/revert_series_lagrange.c +++ b/fmpz_poly/revert_series_lagrange.c @@ -49,7 +49,7 @@ _fmpz_poly_revert_series_lagrange(fmpz * Qinv, const fmpz * Q, slong n) fmpz_zero(Qinv); fmpz_set(Qinv + 1, Q + 1); - _fmpz_poly_inv_series(R, Q + 1, n - 1); + _fmpz_poly_inv_series(R, Q + 1, n - 1, n - 1); _fmpz_vec_set(S, R, n - 1); for (i = 2; i < n; i++) diff --git a/fmpz_poly/revert_series_lagrange_fast.c b/fmpz_poly/revert_series_lagrange_fast.c index 95c165c3ff..6218d1a567 100644 --- a/fmpz_poly/revert_series_lagrange_fast.c +++ b/fmpz_poly/revert_series_lagrange_fast.c @@ -57,7 +57,7 @@ _fmpz_poly_revert_series_lagrange_fast(fmpz * Qinv, const fmpz * Q, slong n) fmpz_zero(Qinv); fmpz_set(Qinv + 1, Q + 1); - _fmpz_poly_inv_series(Ri(1), Q + 1, n - 1); + _fmpz_poly_inv_series(Ri(1), Q + 1, n - 1, n - 1); for (i = 2; i <= m; i++) _fmpz_poly_mullow(Ri(i), Ri(i-1), n - 1, Ri(1), n - 1, n - 1); for (i = 2; i < m; i++) diff --git a/fmpz_poly/test/t-inv_series.c b/fmpz_poly/test/t-inv_series.c new file mode 100644 index 0000000000..33d3be33d5 --- /dev/null +++ b/fmpz_poly/test/t-inv_series.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + FLINT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with FLINT; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2009 William Hart + +******************************************************************************/ + +#include +#include +#include +#include "flint.h" +#include "fmpz.h" +#include "fmpz_poly.h" +#include "ulong_extras.h" + +int +main(void) +{ + int i, result; + FLINT_TEST_INIT(state); + + flint_printf("inv_series...."); + fflush(stdout); + + /* Check Q^{-1} * Q is congruent 1 mod t^n */ + for (i = 0; i < 100 * flint_test_multiplier(); i++) + { + fmpz_poly_t a, b, c, one; + slong n = n_randint(state, 80) + 1; + + fmpz_poly_init(a); + fmpz_poly_init(b); + fmpz_poly_init(c); + fmpz_poly_init(one); + + fmpz_poly_randtest_not_zero(a, state, n_randint(state, 80) + 1, 100); + fmpz_poly_set_coeff_si(a, 0, n_randint(state, 2) ? 1 : -1); + + fmpz_poly_set_ui(one, 1); + + fmpz_poly_inv_series(b, a, n); + fmpz_poly_mullow(c, a, b, n); + + result = (fmpz_poly_equal(c, one)); + if (!result) + { + flint_printf("FAIL:\n"); + flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n\n"); + flint_printf("b = "), fmpz_poly_print(b), flint_printf("\n\n"); + flint_printf("c = "), fmpz_poly_print(c), flint_printf("\n\n"); + abort(); + } + + fmpz_poly_clear(a); + fmpz_poly_clear(b); + fmpz_poly_clear(c); + fmpz_poly_clear(one); + } + + /* Verify bug fix for the case Q = -1 mod (x) */ + { + fmpz_poly_t a, b, c, one; + slong n = 1; + + fmpz_poly_init(a); + fmpz_poly_init(b); + fmpz_poly_init(c); + fmpz_poly_init(one); + + fmpz_poly_set_si(a, -1); + fmpz_poly_set_ui(one, 1); + + fmpz_poly_inv_series(b, a, n); + fmpz_poly_mullow(c, a, b, n); + + result = (fmpz_poly_equal(c, one)); + if (!result) + { + flint_printf("FAIL:\n"); + flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n\n"); + flint_printf("b = "), fmpz_poly_print(b), flint_printf("\n\n"); + flint_printf("c = "), fmpz_poly_print(c), flint_printf("\n\n"); + abort(); + } + + fmpz_poly_clear(a); + fmpz_poly_clear(b); + fmpz_poly_clear(c); + fmpz_poly_clear(one); + } + + FLINT_TEST_CLEANUP(state); + + flint_printf("PASS\n"); + return 0; +} diff --git a/fmpz_poly/test/t-inv_series_basecase.c b/fmpz_poly/test/t-inv_series_basecase.c new file mode 100644 index 0000000000..d059f7e263 --- /dev/null +++ b/fmpz_poly/test/t-inv_series_basecase.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + FLINT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with FLINT; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2009 William Hart + +******************************************************************************/ + +#include +#include +#include +#include "flint.h" +#include "fmpz.h" +#include "fmpz_poly.h" +#include "ulong_extras.h" + +int +main(void) +{ + int i, result; + FLINT_TEST_INIT(state); + + flint_printf("inv_series_basecase...."); + fflush(stdout); + + /* Check Q^{-1} * Q is congruent 1 mod t^n */ + for (i = 0; i < 100 * flint_test_multiplier(); i++) + { + fmpz_poly_t a, b, c, one; + slong n = n_randint(state, 80) + 1; + + fmpz_poly_init(a); + fmpz_poly_init(b); + fmpz_poly_init(c); + fmpz_poly_init(one); + + fmpz_poly_randtest_not_zero(a, state, n_randint(state, 80) + 1, 100); + fmpz_poly_set_coeff_si(a, 0, n_randint(state, 2) ? 1 : -1); + + fmpz_poly_set_ui(one, 1); + + fmpz_poly_inv_series_basecase(b, a, n); + fmpz_poly_mullow(c, a, b, n); + + result = (fmpz_poly_equal(c, one)); + if (!result) + { + flint_printf("FAIL:\n"); + flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n\n"); + flint_printf("b = "), fmpz_poly_print(b), flint_printf("\n\n"); + flint_printf("c = "), fmpz_poly_print(c), flint_printf("\n\n"); + abort(); + } + + fmpz_poly_clear(a); + fmpz_poly_clear(b); + fmpz_poly_clear(c); + fmpz_poly_clear(one); + } + + /* Verify bug fix for the case Q = -1 mod (x) */ + { + fmpz_poly_t a, b, c, one; + slong n = 1; + + fmpz_poly_init(a); + fmpz_poly_init(b); + fmpz_poly_init(c); + fmpz_poly_init(one); + + fmpz_poly_set_si(a, -1); + fmpz_poly_set_ui(one, 1); + + fmpz_poly_inv_series_basecase(b, a, n); + fmpz_poly_mullow(c, a, b, n); + + result = (fmpz_poly_equal(c, one)); + if (!result) + { + flint_printf("FAIL:\n"); + flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n\n"); + flint_printf("b = "), fmpz_poly_print(b), flint_printf("\n\n"); + flint_printf("c = "), fmpz_poly_print(c), flint_printf("\n\n"); + abort(); + } + + fmpz_poly_clear(a); + fmpz_poly_clear(b); + fmpz_poly_clear(c); + fmpz_poly_clear(one); + } + + FLINT_TEST_CLEANUP(state); + + flint_printf("PASS\n"); + return 0; +} diff --git a/fmpz_poly/test/t-inv_series_newton.c b/fmpz_poly/test/t-inv_series_newton.c index c21a7b92c0..a1d086019b 100644 --- a/fmpz_poly/test/t-inv_series_newton.c +++ b/fmpz_poly/test/t-inv_series_newton.c @@ -40,8 +40,6 @@ main(void) flint_printf("inv_series_newton...."); fflush(stdout); - - /* Check Q^{-1} * Q is congruent 1 mod t^n */ for (i = 0; i < 100 * flint_test_multiplier(); i++) { diff --git a/fmpz_poly/theta_qexp.c b/fmpz_poly/theta_qexp.c index 33b915f83c..cb8eecd5f3 100644 --- a/fmpz_poly/theta_qexp.c +++ b/fmpz_poly/theta_qexp.c @@ -64,7 +64,7 @@ _fmpz_poly_theta_qexp(fmpz * f, slong k, slong n) { fmpz * t = _fmpz_vec_init(n); _fmpz_poly_theta_qexp(t, -k, n); - _fmpz_poly_inv_series(f, t, n); + _fmpz_poly_inv_series(f, t, n, n); _fmpz_vec_clear(t, n); return; }