diff --git a/doc/source/fmpcb_poly.rst b/doc/source/fmpcb_poly.rst index 1440b0fa..f0c5ef36 100644 --- a/doc/source/fmpcb_poly.rst +++ b/doc/source/fmpcb_poly.rst @@ -304,6 +304,74 @@ Arithmetic Divides `A` by the polynomial `x - c`, computing the quotient `Q` as well as the remainder `R = f(c)`. +Composition +------------------------------------------------------------------------------- + +.. function:: void _fmpcb_poly_compose_horner(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, fmpcb_srcptr poly2, long len2, long prec) + +.. function:: void fmpcb_poly_compose_horner(fmpcb_poly_t res, const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec) + +.. function:: void _fmpcb_poly_compose_divconquer(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, fmpcb_srcptr poly2, long len2, long prec) + +.. function:: void fmpcb_poly_compose_divconquer(fmpcb_poly_t res, const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec) + +.. function:: void _fmpcb_poly_compose(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, fmpcb_srcptr poly2, long len2, long prec) + +.. function:: void fmpcb_poly_compose(fmpcb_poly_t res, const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec) + + Sets *res* to the composition `h(x) = f(g(x))` where `f` is given by + *poly1* and `g` is given by *poly2*, respectively using Horner's rule, + divide-and-conquer, and an automatic choice between the two algorithms. + The underscore methods do not support aliasing of the output + with either input polynomial. + +.. function:: void _fmpcb_poly_compose_series_horner(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, fmpcb_srcptr poly2, long len2, long n, long prec) + +.. function:: void fmpcb_poly_compose_series_horner(fmpcb_poly_t res, const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long n, long prec) + +.. function:: void _fmpcb_poly_compose_series_brent_kung(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, fmpcb_srcptr poly2, long len2, long n, long prec) + +.. function:: void fmpcb_poly_compose_series_brent_kung(fmpcb_poly_t res, const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long n, long prec) + +.. function:: void _fmpcb_poly_compose_series(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, fmpcb_srcptr poly2, long len2, long n, long prec) + +.. function:: void fmpcb_poly_compose_series(fmpcb_poly_t res, const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long n, long prec) + + Sets *res* to the power series composition `h(x) = f(g(x))` truncated + to order `O(x^n)` where `f` is given by *poly1* and `g` is given by *poly2*, + respectively using Horner's rule, the Brent-Kung baby step-giant step + algorithm, and an automatic choice between the two algorithms. + We require that the constant term in `g(x)` is exactly zero. + The underscore methods do not support aliasing of the output + with either input polynomial. + + +.. function:: void _fmpcb_poly_revert_series_lagrange(fmpcb_ptr h, fmpcb_srcptr f, long n, long prec) + +.. function:: void fmpcb_poly_revert_series_lagrange(fmpcb_poly_t h, const fmpcb_poly_t f, long n, long prec) + +.. function:: void _fmpcb_poly_revert_series_newton(fmpcb_ptr h, fmpcb_srcptr f, long n, long prec) + +.. function:: void fmpcb_poly_revert_series_newton(fmpcb_poly_t h, const fmpcb_poly_t f, long n, long prec) + +.. function:: void _fmpcb_poly_revert_series_lagrange_fast(fmpcb_ptr h, fmpcb_srcptr f, long n, long prec) + +.. function:: void fmpcb_poly_revert_series_lagrange_fast(fmpcb_poly_t h, const fmpcb_poly_t f, long n, long prec) + +.. function:: void _fmpcb_poly_revert_series(fmpcb_ptr h, fmpcb_srcptr f, long n, long prec) + +.. function:: void fmpcb_poly_revert_series(fmpcb_poly_t h, const fmpcb_poly_t f, long n, long prec) + + Sets `h` to the power series reversion of `f`, i.e. the expansion + of the compositional inverse function `f^{-1}(x)`, + truncated to order `O(x^n)`, using respectively + Lagrange inversion, Newton iteration, fast Lagrange inversion, + and a default algorithm choice. + + We require that the constant term in `f` is exactly zero and that the + linear term is nonzero. The underscore methods assume that `f` is zero-padded to length `n` + and do not support aliasing. + Evaluation ------------------------------------------------------------------------------- diff --git a/fmpcb_poly.h b/fmpcb_poly.h index 85a501d3..bb275863 100644 --- a/fmpcb_poly.h +++ b/fmpcb_poly.h @@ -283,6 +283,66 @@ void fmpcb_poly_divrem(fmpcb_poly_t Q, fmpcb_poly_t R, void _fmpcb_poly_div_root(fmpcb_ptr Q, fmpcb_t R, fmpcb_srcptr A, long len, const fmpcb_t c, long prec); +/* Composition */ + +void _fmpcb_poly_compose(fmpcb_ptr res, + fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long prec); + +void fmpcb_poly_compose(fmpcb_poly_t res, + const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec); + +void _fmpcb_poly_compose_horner(fmpcb_ptr res, + fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long prec); + +void fmpcb_poly_compose_horner(fmpcb_poly_t res, + const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec); + +void _fmpcb_poly_compose_divconquer(fmpcb_ptr res, + fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long prec); + +void fmpcb_poly_compose_divconquer(fmpcb_poly_t res, + const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec); + +void _fmpcb_poly_compose_series_horner(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long n, long prec); + +void fmpcb_poly_compose_series_horner(fmpcb_poly_t res, + const fmpcb_poly_t poly1, + const fmpcb_poly_t poly2, long n, long prec); + +void _fmpcb_poly_compose_series_brent_kung(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long n, long prec); + +void fmpcb_poly_compose_series_brent_kung(fmpcb_poly_t res, + const fmpcb_poly_t poly1, + const fmpcb_poly_t poly2, long n, long prec); + +void _fmpcb_poly_compose_series(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long n, long prec); + +void fmpcb_poly_compose_series(fmpcb_poly_t res, + const fmpcb_poly_t poly1, + const fmpcb_poly_t poly2, long n, long prec); + +/* Reversion */ + +void _fmpcb_poly_revert_series_lagrange(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec); +void fmpcb_poly_revert_series_lagrange(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec); + +void _fmpcb_poly_revert_series_newton(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec); +void fmpcb_poly_revert_series_newton(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec); + +void _fmpcb_poly_revert_series_lagrange_fast(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec); +void fmpcb_poly_revert_series_lagrange_fast(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec); + +void _fmpcb_poly_revert_series(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec); +void fmpcb_poly_revert_series(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec); + + + void _fmpcb_poly_evaluate_vec_fast_precomp(fmpcb_ptr vs, fmpcb_srcptr poly, long plen, fmpcb_ptr * tree, long len, long prec); diff --git a/fmpcb_poly/compose.c b/fmpcb_poly/compose.c new file mode 100644 index 00000000..cce67557 --- /dev/null +++ b/fmpcb_poly/compose.c @@ -0,0 +1,78 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2010 William Hart + Copyright (C) 2012 Sebastian Pancratz + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_compose(fmpcb_ptr res, + fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long prec) +{ + if (len1 <= 7) + _fmpcb_poly_compose_horner(res, poly1, len1, poly2, len2, prec); + else + _fmpcb_poly_compose_divconquer(res, poly1, len1, poly2, len2, prec); +} + +void fmpcb_poly_compose(fmpcb_poly_t res, + const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec) +{ + const long len1 = poly1->length; + const long len2 = poly2->length; + + if (len1 == 0) + { + fmpcb_poly_zero(res); + } + else if (len1 == 1 || len2 == 0) + { + fmpcb_poly_set_fmpcb(res, poly1->coeffs); + } + else + { + const long lenr = (len1 - 1) * (len2 - 1) + 1; + + if (res != poly1 && res != poly2) + { + fmpcb_poly_fit_length(res, lenr); + _fmpcb_poly_compose(res->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, lenr); + _fmpcb_poly_compose(t->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, prec); + fmpcb_poly_swap(res, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(res, lenr); + _fmpcb_poly_normalise(res); + } +} diff --git a/fmpcb_poly/compose_divconquer.c b/fmpcb_poly/compose_divconquer.c new file mode 100644 index 00000000..cdfdc752 --- /dev/null +++ b/fmpcb_poly/compose_divconquer.c @@ -0,0 +1,196 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2010 William Hart + Copyright (C) 2012 Sebastian Pancratz + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_compose_divconquer(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long prec) +{ + long i, j, k, n; + long *hlen, alloc, powlen; + fmpcb_ptr v, pow, temp; + fmpcb_ptr * h; + + if (len1 == 1) + { + fmpcb_set(res, poly1); + return; + } + if (len2 == 1) + { + _fmpcb_poly_evaluate(res, poly1, len1, poly2, prec); + return; + } + if (len1 == 2) + { + _fmpcb_poly_compose_horner(res, poly1, len1, poly2, len2, prec); + return; + } + + /* Initialisation */ + + hlen = (long *) flint_malloc(((len1 + 1) / 2) * sizeof(long)); + + for (k = 1; (2 << k) < len1; k++) ; + + hlen[0] = hlen[1] = ((1 << k) - 1) * (len2 - 1) + 1; + for (i = k - 1; i > 0; i--) + { + long hi = (len1 + (1 << i) - 1) / (1 << i); + for (n = (hi + 1) / 2; n < hi; n++) + hlen[n] = ((1 << i) - 1) * (len2 - 1) + 1; + } + powlen = (1 << k) * (len2 - 1) + 1; + + alloc = 0; + for (i = 0; i < (len1 + 1) / 2; i++) + alloc += hlen[i]; + + v = _fmpcb_vec_init(alloc + 2 * powlen); + h = (fmpcb_ptr *) flint_malloc(((len1 + 1) / 2) * sizeof(fmpcb_ptr)); + h[0] = v; + for (i = 0; i < (len1 - 1) / 2; i++) + { + h[i + 1] = h[i] + hlen[i]; + hlen[i] = 0; + } + hlen[(len1 - 1) / 2] = 0; + pow = v + alloc; + temp = pow + powlen; + + /* Let's start the actual work */ + + for (i = 0, j = 0; i < len1 / 2; i++, j += 2) + { + if (!fmpcb_is_zero(poly1 + j + 1)) + { + _fmpcb_vec_scalar_mul(h[i], poly2, len2, poly1 + j + 1, prec); + fmpcb_add(h[i], h[i], poly1 + j, prec); + hlen[i] = len2; + } + else if (!fmpcb_is_zero(poly1 + j)) + { + fmpcb_set(h[i], poly1 + j); + hlen[i] = 1; + } + } + if ((len1 & 1L)) + { + if (!fmpcb_is_zero(poly1 + j)) + { + fmpcb_set(h[i], poly1 + j); + hlen[i] = 1; + } + } + + _fmpcb_poly_mul(pow, poly2, len2, poly2, len2, prec); + powlen = 2 * len2 - 1; + + for (n = (len1 + 1) / 2; n > 2; n = (n + 1) / 2) + { + if (hlen[1] > 0) + { + long templen = powlen + hlen[1] - 1; + _fmpcb_poly_mul(temp, pow, powlen, h[1], hlen[1], prec); + _fmpcb_poly_add(h[0], temp, templen, h[0], hlen[0], prec); + hlen[0] = FLINT_MAX(hlen[0], templen); + } + + for (i = 1; i < n / 2; i++) + { + if (hlen[2*i + 1] > 0) + { + _fmpcb_poly_mul(h[i], pow, powlen, h[2*i + 1], hlen[2*i + 1], prec); + hlen[i] = hlen[2*i + 1] + powlen - 1; + } else + hlen[i] = 0; + _fmpcb_poly_add(h[i], h[i], hlen[i], h[2*i], hlen[2*i], prec); + hlen[i] = FLINT_MAX(hlen[i], hlen[2*i]); + } + if ((n & 1L)) + { + _fmpcb_vec_set(h[i], h[2*i], hlen[2*i]); + hlen[i] = hlen[2*i]; + } + + _fmpcb_poly_mul(temp, pow, powlen, pow, powlen, prec); + powlen += powlen - 1; + { + fmpcb_ptr t = pow; + pow = temp; + temp = t; + } + } + + _fmpcb_poly_mul(res, pow, powlen, h[1], hlen[1], prec); + _fmpcb_vec_add(res, res, h[0], hlen[0], prec); + + _fmpcb_vec_clear(v, alloc + 2 * powlen); + flint_free(h); + flint_free(hlen); +} + +void +fmpcb_poly_compose_divconquer(fmpcb_poly_t res, + const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec) +{ + const long len1 = poly1->length; + const long len2 = poly2->length; + + if (len1 == 0) + { + fmpcb_poly_zero(res); + } + else if (len1 == 1 || len2 == 0) + { + fmpcb_poly_set_fmpcb(res, poly1->coeffs); + } + else + { + const long lenr = (len1 - 1) * (len2 - 1) + 1; + + if (res != poly1 && res != poly2) + { + fmpcb_poly_fit_length(res, lenr); + _fmpcb_poly_compose_divconquer(res->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, lenr); + _fmpcb_poly_compose_divconquer(t->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, prec); + fmpcb_poly_swap(res, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(res, lenr); + _fmpcb_poly_normalise(res); + } +} diff --git a/fmpcb_poly/compose_horner.c b/fmpcb_poly/compose_horner.c new file mode 100644 index 00000000..1a6f00ed --- /dev/null +++ b/fmpcb_poly/compose_horner.c @@ -0,0 +1,125 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2010 William Hart + Copyright (C) 2012 Sebastian Pancratz + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_compose_horner(fmpcb_ptr res, + fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long prec) +{ + if (len1 == 1) + { + fmpcb_set(res, poly1); + } + else if (len2 == 1) + { + _fmpcb_poly_evaluate(res, poly1, len1, poly2, prec); + } + else if (len1 == 2) + { + _fmpcb_vec_scalar_mul(res, poly2, len2, poly1 + 1, prec); + fmpcb_add(res, res, poly1, prec); + } + else + { + const long alloc = (len1 - 1) * (len2 - 1) + 1; + long i = len1 - 1, lenr = len2; + fmpcb_ptr t, t1, t2; + t = _fmpcb_vec_init(alloc); + + if (len1 % 2 == 0) + { + t1 = res; + t2 = t; + } + else + { + t1 = t; + t2 = res; + } + + /* Perform the first two steps as one, + "res = a(m) * poly2 + a(m-1)". */ + { + _fmpcb_vec_scalar_mul(t1, poly2, len2, poly1 + i, prec); + i--; + fmpcb_add(t1 + 0, t1 + 0, poly1 + i, prec); + } + while (i--) + { + _fmpcb_poly_mul(t2, t1, lenr, poly2, len2, prec); + lenr += len2 - 1; + { + void *t_ = t1; + t1 = t2; + t2 = t_; + } + fmpcb_add(t1 + 0, t1 + 0, poly1 + i, prec); + } + _fmpcb_vec_clear(t, alloc); + } +} + +void fmpcb_poly_compose_horner(fmpcb_poly_t res, + const fmpcb_poly_t poly1, const fmpcb_poly_t poly2, long prec) +{ + const long len1 = poly1->length; + const long len2 = poly2->length; + + if (len1 == 0) + { + fmpcb_poly_zero(res); + } + else if (len1 == 1 || len2 == 0) + { + fmpcb_poly_set_fmpcb(res, poly1->coeffs); + } + else + { + const long lenr = (len1 - 1) * (len2 - 1) + 1; + + if (res != poly1 && res != poly2) + { + fmpcb_poly_fit_length(res, lenr); + _fmpcb_poly_compose_horner(res->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, lenr); + _fmpcb_poly_compose_horner(t->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, prec); + fmpcb_poly_swap(res, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(res, lenr); + _fmpcb_poly_normalise(res); + } +} diff --git a/fmpcb_poly/compose_series.c b/fmpcb_poly/compose_series.c new file mode 100644 index 00000000..bd4dd37e --- /dev/null +++ b/fmpcb_poly/compose_series.c @@ -0,0 +1,112 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_compose_series(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long n, long prec) +{ + if (len2 == 2) /* fast linear case (TODO: all monomials) */ + { + long i; + fmpcb_t t; + + fmpcb_init(t); + + fmpcb_set(t, poly2 + 1); + fmpcb_set_round(res, poly1, prec); + + for (i = 1; i < n; i++) + { + fmpcb_mul(res + i, poly1 + i, t, prec); + if (i + 1 < n) + fmpcb_mul(t, t, poly2 + 1, prec); + } + + fmpcb_clear(t); + } + else if (len1 < 6 || n < 6) + { + _fmpcb_poly_compose_series_horner(res, poly1, len1, poly2, len2, n, prec); + } + else + { + _fmpcb_poly_compose_series_brent_kung(res, poly1, len1, poly2, len2, n, prec); + } +} + +void +fmpcb_poly_compose_series(fmpcb_poly_t res, + const fmpcb_poly_t poly1, + const fmpcb_poly_t poly2, long n, long prec) +{ + long len1 = poly1->length; + long len2 = poly2->length; + long lenr; + + if (len2 != 0 && !fmpcb_is_zero(poly2->coeffs)) + { + printf("exception: compose_series: inner " + "polynomial must have zero constant term\n"); + abort(); + } + + if (len1 == 0 || n == 0) + { + fmpcb_poly_zero(res); + return; + } + + if (len2 == 0 || len1 == 1) + { + fmpcb_poly_set_fmpcb(res, poly1->coeffs); + return; + } + + lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n); + len1 = FLINT_MIN(len1, lenr); + len2 = FLINT_MIN(len2, lenr); + + if ((res != poly1) && (res != poly2)) + { + fmpcb_poly_fit_length(res, lenr); + _fmpcb_poly_compose_series(res->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, lenr, prec); + _fmpcb_poly_set_length(res, lenr); + _fmpcb_poly_normalise(res); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, lenr); + _fmpcb_poly_compose_series(t->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, lenr, prec); + _fmpcb_poly_set_length(t, lenr); + _fmpcb_poly_normalise(t); + fmpcb_poly_swap(res, t); + fmpcb_poly_clear(t); + } +} diff --git a/fmpcb_poly/compose_series_brent_kung.c b/fmpcb_poly/compose_series_brent_kung.c new file mode 100644 index 00000000..052fbf94 --- /dev/null +++ b/fmpcb_poly/compose_series_brent_kung.c @@ -0,0 +1,135 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" +#include "fmpcb_mat.h" + +void +_fmpcb_poly_compose_series_brent_kung(fmpcb_ptr res, + fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long n, long prec) +{ + fmpcb_mat_t A, B, C; + fmpcb_ptr t, h; + long i, m; + + if (n == 1) + { + fmpcb_set(res, poly1); + return; + } + + m = n_sqrt(n) + 1; + + fmpcb_mat_init(A, m, n); + fmpcb_mat_init(B, m, m); + fmpcb_mat_init(C, m, n); + + h = _fmpcb_vec_init(n); + t = _fmpcb_vec_init(n); + + /* Set rows of B to the segments of poly1 */ + for (i = 0; i < len1 / m; i++) + _fmpcb_vec_set(B->rows[i], poly1 + i*m, m); + _fmpcb_vec_set(B->rows[i], poly1 + i*m, len1 % m); + + /* Set rows of A to powers of poly2 */ + fmpcb_set_ui(A->rows[0] + 0, 1UL); + _fmpcb_vec_set(A->rows[1], poly2, len2); + for (i = 2; i < m; i++) + _fmpcb_poly_mullow(A->rows[i], A->rows[(i + 1) / 2], n, A->rows[i / 2], n, n, prec); + + fmpcb_mat_mul(C, B, A, prec); + + /* Evaluate block composition using the Horner scheme */ + _fmpcb_vec_set(res, C->rows[m - 1], n); + _fmpcb_poly_mullow(h, A->rows[m - 1], n, poly2, len2, n, prec); + + for (i = m - 2; i >= 0; i--) + { + _fmpcb_poly_mullow(t, res, n, h, n, n, prec); + _fmpcb_poly_add(res, t, n, C->rows[i], n, prec); + } + + _fmpcb_vec_clear(h, n); + _fmpcb_vec_clear(t, n); + + fmpcb_mat_clear(A); + fmpcb_mat_clear(B); + fmpcb_mat_clear(C); +} + +void +fmpcb_poly_compose_series_brent_kung(fmpcb_poly_t res, + const fmpcb_poly_t poly1, + const fmpcb_poly_t poly2, long n, long prec) +{ + long len1 = poly1->length; + long len2 = poly2->length; + long lenr; + + if (len2 != 0 && !fmpcb_is_zero(poly2->coeffs)) + { + printf("exception: compose_series: inner " + "polynomial must have zero constant term\n"); + abort(); + } + + if (len1 == 0 || n == 0) + { + fmpcb_poly_zero(res); + return; + } + + if (len2 == 0 || len1 == 1) + { + fmpcb_poly_set_fmpcb(res, poly1->coeffs); + return; + } + + lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n); + len1 = FLINT_MIN(len1, lenr); + len2 = FLINT_MIN(len2, lenr); + + if ((res != poly1) && (res != poly2)) + { + fmpcb_poly_fit_length(res, lenr); + _fmpcb_poly_compose_series_brent_kung(res->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, lenr, prec); + _fmpcb_poly_set_length(res, lenr); + _fmpcb_poly_normalise(res); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, lenr); + _fmpcb_poly_compose_series_brent_kung(t->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, lenr, prec); + _fmpcb_poly_set_length(t, lenr); + _fmpcb_poly_normalise(t); + fmpcb_poly_swap(res, t); + fmpcb_poly_clear(t); + } +} diff --git a/fmpcb_poly/compose_series_horner.c b/fmpcb_poly/compose_series_horner.c new file mode 100644 index 00000000..87de140b --- /dev/null +++ b/fmpcb_poly/compose_series_horner.c @@ -0,0 +1,120 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_compose_series_horner(fmpcb_ptr res, fmpcb_srcptr poly1, long len1, + fmpcb_srcptr poly2, long len2, long n, long prec) +{ + if (n == 1) + { + fmpcb_set(res, poly1); + } + else + { + long i = len1 - 1; + long lenr; + + fmpcb_ptr t = _fmpcb_vec_init(n); + + lenr = len2; + _fmpcb_vec_scalar_mul(res, poly2, len2, poly1 + i, prec); + i--; + fmpcb_add(res, res, poly1 + i, prec); + + while (i > 0) + { + i--; + if (lenr + len2 - 1 < n) + { + _fmpcb_poly_mul(t, res, lenr, poly2, len2, prec); + lenr = lenr + len2 - 1; + } + else + { + _fmpcb_poly_mullow(t, res, lenr, poly2, len2, n, prec); + lenr = n; + } + _fmpcb_poly_add(res, t, lenr, poly1 + i, 1, prec); + } + + _fmpcb_vec_zero(res + lenr, n - lenr); + _fmpcb_vec_clear(t, n); + } +} + +void +fmpcb_poly_compose_series_horner(fmpcb_poly_t res, + const fmpcb_poly_t poly1, + const fmpcb_poly_t poly2, long n, long prec) +{ + long len1 = poly1->length; + long len2 = poly2->length; + long lenr; + + if (len2 != 0 && !fmpcb_is_zero(poly2->coeffs)) + { + printf("exception: compose_series: inner " + "polynomial must have zero constant term\n"); + abort(); + } + + if (len1 == 0 || n == 0) + { + fmpcb_poly_zero(res); + return; + } + + if (len2 == 0 || len1 == 1) + { + fmpcb_poly_set_fmpcb(res, poly1->coeffs); + return; + } + + lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n); + len1 = FLINT_MIN(len1, lenr); + len2 = FLINT_MIN(len2, lenr); + + if ((res != poly1) && (res != poly2)) + { + fmpcb_poly_fit_length(res, lenr); + _fmpcb_poly_compose_series_horner(res->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, lenr, prec); + _fmpcb_poly_set_length(res, lenr); + _fmpcb_poly_normalise(res); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, lenr); + _fmpcb_poly_compose_series_horner(t->coeffs, poly1->coeffs, len1, + poly2->coeffs, len2, lenr, prec); + _fmpcb_poly_set_length(t, lenr); + _fmpcb_poly_normalise(t); + fmpcb_poly_swap(res, t); + fmpcb_poly_clear(t); + } +} diff --git a/fmpcb_poly/revert_series.c b/fmpcb_poly/revert_series.c new file mode 100644 index 00000000..073c68ad --- /dev/null +++ b/fmpcb_poly/revert_series.c @@ -0,0 +1,91 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_revert_series(fmpcb_ptr Qinv, + fmpcb_srcptr Q, long n, long prec) +{ + _fmpcb_poly_revert_series_lagrange_fast(Qinv, Q, n, prec); +} + +void +fmpcb_poly_revert_series(fmpcb_poly_t Qinv, + const fmpcb_poly_t Q, long n, long prec) +{ + fmpcb_ptr Qcopy; + int Qalloc; + long Qlen = Q->length; + + if (Q->length < 2 || !fmpcb_is_zero(Q->coeffs) + || fmpcb_contains_zero(Q->coeffs + 1)) + { + printf("Exception (fmpcb_poly_revert_series). Input must \n" + "have zero constant term and nonzero coefficient of x^1.\n"); + abort(); + } + + if (n < 2) + { + fmpcb_poly_zero(Qinv); + return; + } + + if (Qlen >= n) + { + Qcopy = Q->coeffs; + Qalloc = 0; + } + else + { + long i; + Qcopy = _fmpcb_vec_init(n); + for (i = 0; i < Qlen; i++) + Qcopy[i] = Q->coeffs[i]; + Qalloc = 1; + } + + if (Qinv != Q) + { + fmpcb_poly_fit_length(Qinv, n); + _fmpcb_poly_revert_series(Qinv->coeffs, Qcopy, n, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, n); + _fmpcb_poly_revert_series(t->coeffs, Qcopy, n, prec); + fmpcb_poly_swap(Qinv, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(Qinv, n); + _fmpcb_poly_normalise(Qinv); + + if (Qalloc) + flint_free(Qcopy); +} + diff --git a/fmpcb_poly/revert_series_lagrange.c b/fmpcb_poly/revert_series_lagrange.c new file mode 100644 index 00000000..f58db6bc --- /dev/null +++ b/fmpcb_poly/revert_series_lagrange.c @@ -0,0 +1,122 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +void +_fmpcb_poly_revert_series_lagrange(fmpcb_ptr Qinv, + fmpcb_srcptr Q, long n, long prec) +{ + long i; + fmpcb_ptr R, S, T, tmp; + + if (n <= 2) + { + if (n >= 1) + fmpcb_zero(Qinv); + if (n == 2) + fmpcb_inv(Qinv + 1, Q + 1, prec); + return; + } + + R = _fmpcb_vec_init(n - 1); + S = _fmpcb_vec_init(n - 1); + T = _fmpcb_vec_init(n - 1); + + fmpcb_zero(Qinv); + fmpcb_inv(Qinv + 1, Q + 1, prec); + + _fmpcb_poly_inv_series(R, Q + 1, n - 1, n - 1, prec); + _fmpcb_vec_set(S, R, n - 1); + + for (i = 2; i < n; i++) + { + _fmpcb_poly_mullow(T, S, n - 1, R, n - 1, n - 1, prec); + fmpcb_div_ui(Qinv + i, T + i - 1, i, prec); + tmp = S; S = T; T = tmp; + } + + _fmpcb_vec_clear(R, n - 1); + _fmpcb_vec_clear(S, n - 1); + _fmpcb_vec_clear(T, n - 1); +} + +void +fmpcb_poly_revert_series_lagrange(fmpcb_poly_t Qinv, + const fmpcb_poly_t Q, long n, long prec) +{ + fmpcb_ptr Qcopy; + int Qalloc; + long Qlen = Q->length; + + if (Q->length < 2 || !fmpcb_is_zero(Q->coeffs) + || fmpcb_contains_zero(Q->coeffs + 1)) + { + printf("Exception (fmpcb_poly_revert_series_lagrange). Input must \n" + "have zero constant term and nonzero coefficient of x^1.\n"); + abort(); + } + + if (n < 2) + { + fmpcb_poly_zero(Qinv); + return; + } + + if (Qlen >= n) + { + Qcopy = Q->coeffs; + Qalloc = 0; + } + else + { + long i; + Qcopy = _fmpcb_vec_init(n); + for (i = 0; i < Qlen; i++) + Qcopy[i] = Q->coeffs[i]; + Qalloc = 1; + } + + if (Qinv != Q) + { + fmpcb_poly_fit_length(Qinv, n); + _fmpcb_poly_revert_series_lagrange(Qinv->coeffs, Qcopy, n, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, n); + _fmpcb_poly_revert_series_lagrange(t->coeffs, Qcopy, n, prec); + fmpcb_poly_swap(Qinv, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(Qinv, n); + _fmpcb_poly_normalise(Qinv); + + if (Qalloc) + flint_free(Qcopy); +} + diff --git a/fmpcb_poly/revert_series_lagrange_fast.c b/fmpcb_poly/revert_series_lagrange_fast.c new file mode 100644 index 00000000..f3f182d5 --- /dev/null +++ b/fmpcb_poly/revert_series_lagrange_fast.c @@ -0,0 +1,147 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +/* pointer to (x/Q)^i */ +#define Ri(ii) (R + (n-1)*((ii)-1)) + +void +_fmpcb_poly_revert_series_lagrange_fast(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec) +{ + long i, j, k, m; + fmpcb_ptr R, S, T, tmp; + fmpcb_t t; + + if (n <= 2) + { + if (n >= 1) + fmpcb_zero(Qinv); + if (n == 2) + fmpcb_inv(Qinv + 1, Q + 1, prec); + return; + } + + m = n_sqrt(n); + + fmpcb_init(t); + R = _fmpcb_vec_init((n - 1) * m); + S = _fmpcb_vec_init(n - 1); + T = _fmpcb_vec_init(n - 1); + + fmpcb_zero(Qinv); + fmpcb_inv(Qinv + 1, Q + 1, prec); + + _fmpcb_poly_inv_series(Ri(1), Q + 1, n - 1, n - 1, prec); + for (i = 2; i <= m; i++) + _fmpcb_poly_mullow(Ri(i), Ri((i + 1) / 2), n - 1, Ri(i / 2), n - 1, n - 1, prec); + + for (i = 2; i < m; i++) + fmpcb_div_ui(Qinv + i, Ri(i) + i - 1, i, prec); + + _fmpcb_vec_set(S, Ri(m), n - 1); + + for (i = m; i < n; i += m) + { + fmpcb_div_ui(Qinv + i, S + i - 1, i, prec); + + for (j = 1; j < m && i + j < n; j++) + { + fmpcb_mul(t, S + 0, Ri(j) + i + j - 1, prec); + for (k = 1; k <= i + j - 1; k++) + fmpcb_addmul(t, S + k, Ri(j) + i + j - 1 - k, prec); + fmpcb_div_ui(Qinv + i + j, t, i + j, prec); + } + + if (i + 1 < n) + { + _fmpcb_poly_mullow(T, S, n - 1, Ri(m), n - 1, n - 1, prec); + tmp = S; S = T; T = tmp; + } + } + + fmpcb_clear(t); + _fmpcb_vec_clear(R, (n - 1) * m); + _fmpcb_vec_clear(S, n - 1); + _fmpcb_vec_clear(T, n - 1); +} + +void +fmpcb_poly_revert_series_lagrange_fast(fmpcb_poly_t Qinv, + const fmpcb_poly_t Q, long n, long prec) +{ + fmpcb_ptr Qcopy; + int Qalloc; + long Qlen = Q->length; + + if (Q->length < 2 || !fmpcb_is_zero(Q->coeffs) + || fmpcb_contains_zero(Q->coeffs + 1)) + { + printf("Exception (fmpcb_poly_revert_series_lagrange_fast). Input \n" + "must have zero constant term and nonzero coefficient of x^1.\n"); + abort(); + } + + if (n < 2) + { + fmpcb_poly_zero(Qinv); + return; + } + + if (Qlen >= n) + { + Qcopy = Q->coeffs; + Qalloc = 0; + } + else + { + long i; + Qcopy = _fmpcb_vec_init(n); + for (i = 0; i < Qlen; i++) + Qcopy[i] = Q->coeffs[i]; + Qalloc = 1; + } + + if (Qinv != Q) + { + fmpcb_poly_fit_length(Qinv, n); + _fmpcb_poly_revert_series_lagrange_fast(Qinv->coeffs, Qcopy, n, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, n); + _fmpcb_poly_revert_series_lagrange_fast(t->coeffs, Qcopy, n, prec); + fmpcb_poly_swap(Qinv, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(Qinv, n); + _fmpcb_poly_normalise(Qinv); + + if (Qalloc) + flint_free(Qcopy); +} + diff --git a/fmpcb_poly/revert_series_newton.c b/fmpcb_poly/revert_series_newton.c new file mode 100644 index 00000000..d49e6560 --- /dev/null +++ b/fmpcb_poly/revert_series_newton.c @@ -0,0 +1,131 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +#define CUTOFF 5 + +void +_fmpcb_poly_revert_series_newton(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec) +{ + long i, k, a[FLINT_BITS]; + fmpcb_ptr T, U, V; + + if (n <= 2) + { + if (n >= 1) + fmpcb_zero(Qinv); + if (n == 2) + fmpcb_inv(Qinv + 1, Q + 1, prec); + return; + } + + T = _fmpcb_vec_init(n); + U = _fmpcb_vec_init(n); + V = _fmpcb_vec_init(n); + + k = n; + for (i = 1; (1L << i) < k; i++); + a[i = 0] = k; + while (k >= CUTOFF) + a[++i] = (k = (k + 1) / 2); + + _fmpcb_poly_revert_series_lagrange(Qinv, Q, k, prec); + _fmpcb_vec_zero(Qinv + k, n - k); + + for (i--; i >= 0; i--) + { + k = a[i]; + _fmpcb_poly_compose_series(T, Q, k, Qinv, k, k, prec); + _fmpcb_poly_derivative(U, T, k, prec); fmpcb_zero(U + k - 1); + fmpcb_zero(T + 1); + _fmpcb_poly_div_series(V, T, k, U, k, k, prec); + _fmpcb_poly_derivative(T, Qinv, k, prec); + _fmpcb_poly_mullow(U, V, k, T, k, k, prec); + _fmpcb_vec_sub(Qinv, Qinv, U, k, prec); + } + + _fmpcb_vec_clear(T, n); + _fmpcb_vec_clear(U, n); + _fmpcb_vec_clear(V, n); +} + +void +fmpcb_poly_revert_series_newton(fmpcb_poly_t Qinv, + const fmpcb_poly_t Q, long n, long prec) +{ + fmpcb_ptr Qcopy; + int Qalloc; + long Qlen = Q->length; + + if (Q->length < 2 || !fmpcb_is_zero(Q->coeffs) + || fmpcb_contains_zero(Q->coeffs + 1)) + { + printf("Exception (fmpcb_poly_revert_series_newton). Input must \n" + "have zero constant term and nonzero coefficient of x^1.\n"); + abort(); + } + + if (n < 2) + { + fmpcb_poly_zero(Qinv); + return; + } + + if (Qlen >= n) + { + Qcopy = Q->coeffs; + Qalloc = 0; + } + else + { + long i; + Qcopy = _fmpcb_vec_init(n); + for (i = 0; i < Qlen; i++) + Qcopy[i] = Q->coeffs[i]; + Qalloc = 1; + } + + if (Qinv != Q) + { + fmpcb_poly_fit_length(Qinv, n); + _fmpcb_poly_revert_series_newton(Qinv->coeffs, Qcopy, n, prec); + } + else + { + fmpcb_poly_t t; + fmpcb_poly_init2(t, n); + _fmpcb_poly_revert_series_newton(t->coeffs, Qcopy, n, prec); + fmpcb_poly_swap(Qinv, t); + fmpcb_poly_clear(t); + } + + _fmpcb_poly_set_length(Qinv, n); + _fmpcb_poly_normalise(Qinv); + + if (Qalloc) + flint_free(Qcopy); +} + diff --git a/fmpcb_poly/test/t-compose.c b/fmpcb_poly/test/t-compose.c new file mode 100644 index 00000000..cfe1c073 --- /dev/null +++ b/fmpcb_poly/test/t-compose.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("compose...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + long qbits1, qbits2, rbits1, rbits2, rbits3; + fmpq_poly_t A, B, C; + fmpcb_poly_t a, b, c, d; + + qbits1 = 2 + n_randint(state, 200); + qbits2 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + + fmpq_poly_init(A); + fmpq_poly_init(B); + fmpq_poly_init(C); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + fmpcb_poly_init(d); + + fmpq_poly_randtest(A, state, 1 + n_randint(state, 20), qbits1); + fmpq_poly_randtest(B, state, 1 + n_randint(state, 10), qbits2); + fmpq_poly_compose(C, A, B); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_set_fmpq_poly(b, B, rbits2); + fmpcb_poly_compose(c, a, b, rbits3); + + if (!fmpcb_poly_contains_fmpq_poly(c, C)) + { + printf("FAIL\n\n"); + printf("bits3 = %ld\n", rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + printf("C = "); fmpq_poly_print(C); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(d, a); + fmpcb_poly_compose(d, d, b, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmpcb_poly_set(d, b); + fmpcb_poly_compose(d, a, d, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + fmpq_poly_clear(C); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + fmpcb_poly_clear(d); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-compose_divconquer.c b/fmpcb_poly/test/t-compose_divconquer.c new file mode 100644 index 00000000..68086ff8 --- /dev/null +++ b/fmpcb_poly/test/t-compose_divconquer.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("compose_divconquer...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + long qbits1, qbits2, rbits1, rbits2, rbits3; + fmpq_poly_t A, B, C; + fmpcb_poly_t a, b, c, d; + + qbits1 = 2 + n_randint(state, 200); + qbits2 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + + fmpq_poly_init(A); + fmpq_poly_init(B); + fmpq_poly_init(C); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + fmpcb_poly_init(d); + + fmpq_poly_randtest(A, state, 1 + n_randint(state, 20), qbits1); + fmpq_poly_randtest(B, state, 1 + n_randint(state, 10), qbits2); + fmpq_poly_compose(C, A, B); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_set_fmpq_poly(b, B, rbits2); + fmpcb_poly_compose_divconquer(c, a, b, rbits3); + + if (!fmpcb_poly_contains_fmpq_poly(c, C)) + { + printf("FAIL\n\n"); + printf("bits3 = %ld\n", rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + printf("C = "); fmpq_poly_print(C); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(d, a); + fmpcb_poly_compose_divconquer(d, d, b, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmpcb_poly_set(d, b); + fmpcb_poly_compose_divconquer(d, a, d, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + fmpq_poly_clear(C); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + fmpcb_poly_clear(d); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-compose_horner.c b/fmpcb_poly/test/t-compose_horner.c new file mode 100644 index 00000000..376a57ad --- /dev/null +++ b/fmpcb_poly/test/t-compose_horner.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("compose_horner...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + long qbits1, qbits2, rbits1, rbits2, rbits3; + fmpq_poly_t A, B, C; + fmpcb_poly_t a, b, c, d; + + qbits1 = 2 + n_randint(state, 200); + qbits2 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + + fmpq_poly_init(A); + fmpq_poly_init(B); + fmpq_poly_init(C); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + fmpcb_poly_init(d); + + fmpq_poly_randtest(A, state, 1 + n_randint(state, 20), qbits1); + fmpq_poly_randtest(B, state, 1 + n_randint(state, 10), qbits2); + fmpq_poly_compose(C, A, B); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_set_fmpq_poly(b, B, rbits2); + fmpcb_poly_compose_horner(c, a, b, rbits3); + + if (!fmpcb_poly_contains_fmpq_poly(c, C)) + { + printf("FAIL\n\n"); + printf("bits3 = %ld\n", rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + printf("C = "); fmpq_poly_print(C); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(d, a); + fmpcb_poly_compose_horner(d, d, b, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmpcb_poly_set(d, b); + fmpcb_poly_compose_horner(d, a, d, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + fmpq_poly_clear(C); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + fmpcb_poly_clear(d); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-compose_series.c b/fmpcb_poly/test/t-compose_series.c new file mode 100644 index 00000000..905fb6f2 --- /dev/null +++ b/fmpcb_poly/test/t-compose_series.c @@ -0,0 +1,116 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("compose_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + long qbits1, qbits2, rbits1, rbits2, rbits3, n; + fmpq_poly_t A, B, C; + fmpcb_poly_t a, b, c, d; + + qbits1 = 2 + n_randint(state, 200); + qbits2 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + fmpq_poly_init(C); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + fmpcb_poly_init(d); + + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_randtest(B, state, 1 + n_randint(state, 25), qbits2); + fmpq_poly_set_coeff_ui(B, 0, 0); + fmpq_poly_compose_series(C, A, B, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_set_fmpq_poly(b, B, rbits2); + fmpcb_poly_compose_series(c, a, b, n, rbits3); + + if (!fmpcb_poly_contains_fmpq_poly(c, C)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits3 = %ld\n", n, rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + printf("C = "); fmpq_poly_print(C); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(d, a); + fmpcb_poly_compose_series(d, d, b, n, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmpcb_poly_set(d, b); + fmpcb_poly_compose_series(d, a, d, n, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + fmpq_poly_clear(C); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + fmpcb_poly_clear(d); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-compose_series_brent_kung.c b/fmpcb_poly/test/t-compose_series_brent_kung.c new file mode 100644 index 00000000..b3d3c85d --- /dev/null +++ b/fmpcb_poly/test/t-compose_series_brent_kung.c @@ -0,0 +1,116 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("compose_series_brent_kung...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + long qbits1, qbits2, rbits1, rbits2, rbits3, n; + fmpq_poly_t A, B, C; + fmpcb_poly_t a, b, c, d; + + qbits1 = 2 + n_randint(state, 200); + qbits2 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + fmpq_poly_init(C); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + fmpcb_poly_init(d); + + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_randtest(B, state, 1 + n_randint(state, 25), qbits2); + fmpq_poly_set_coeff_ui(B, 0, 0); + fmpq_poly_compose_series(C, A, B, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_set_fmpq_poly(b, B, rbits2); + fmpcb_poly_compose_series_brent_kung(c, a, b, n, rbits3); + + if (!fmpcb_poly_contains_fmpq_poly(c, C)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits3 = %ld\n", n, rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + printf("C = "); fmpq_poly_print(C); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(d, a); + fmpcb_poly_compose_series_brent_kung(d, d, b, n, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmpcb_poly_set(d, b); + fmpcb_poly_compose_series_brent_kung(d, a, d, n, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + fmpq_poly_clear(C); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + fmpcb_poly_clear(d); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-compose_series_horner.c b/fmpcb_poly/test/t-compose_series_horner.c new file mode 100644 index 00000000..3b259e76 --- /dev/null +++ b/fmpcb_poly/test/t-compose_series_horner.c @@ -0,0 +1,116 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("compose_series_horner...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000; iter++) + { + long qbits1, qbits2, rbits1, rbits2, rbits3, n; + fmpq_poly_t A, B, C; + fmpcb_poly_t a, b, c, d; + + qbits1 = 2 + n_randint(state, 200); + qbits2 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + fmpq_poly_init(C); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + fmpcb_poly_init(d); + + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_randtest(B, state, 1 + n_randint(state, 25), qbits2); + fmpq_poly_set_coeff_ui(B, 0, 0); + fmpq_poly_compose_series(C, A, B, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_set_fmpq_poly(b, B, rbits2); + fmpcb_poly_compose_series_horner(c, a, b, n, rbits3); + + if (!fmpcb_poly_contains_fmpq_poly(c, C)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits3 = %ld\n", n, rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + printf("C = "); fmpq_poly_print(C); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(d, a); + fmpcb_poly_compose_series_horner(d, d, b, n, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmpcb_poly_set(d, b); + fmpcb_poly_compose_series_horner(d, a, d, n, rbits3); + if (!fmpcb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + fmpq_poly_clear(C); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + fmpcb_poly_clear(d); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-revert_series.c b/fmpcb_poly/test/t-revert_series.c new file mode 100644 index 00000000..87d40f02 --- /dev/null +++ b/fmpcb_poly/test/t-revert_series.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("revert_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long qbits1, rbits1, rbits2, n; + fmpq_poly_t A, B; + fmpcb_poly_t a, b, c; + + qbits1 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + + do { + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_set_coeff_ui(A, 0, 0); + } while (A->length < 2 || fmpz_is_zero(A->coeffs + 1)); + + fmpq_poly_revert_series(B, A, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_revert_series(b, a, n, rbits2); + + if (!fmpcb_poly_contains_fmpq_poly(b, B)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits2 = %ld\n", n, rbits2); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(c, a); + fmpcb_poly_revert_series(c, c, n, rbits2); + if (!fmpcb_poly_equal(c, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-revert_series_lagrange.c b/fmpcb_poly/test/t-revert_series_lagrange.c new file mode 100644 index 00000000..fb4e4553 --- /dev/null +++ b/fmpcb_poly/test/t-revert_series_lagrange.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("revert_series_lagrange...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long qbits1, rbits1, rbits2, n; + fmpq_poly_t A, B; + fmpcb_poly_t a, b, c; + + qbits1 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + + do { + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_set_coeff_ui(A, 0, 0); + } while (A->length < 2 || fmpz_is_zero(A->coeffs + 1)); + + fmpq_poly_revert_series(B, A, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_revert_series_lagrange(b, a, n, rbits2); + + if (!fmpcb_poly_contains_fmpq_poly(b, B)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits2 = %ld\n", n, rbits2); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(c, a); + fmpcb_poly_revert_series_lagrange(c, c, n, rbits2); + if (!fmpcb_poly_equal(c, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-revert_series_lagrange_fast.c b/fmpcb_poly/test/t-revert_series_lagrange_fast.c new file mode 100644 index 00000000..40c2a920 --- /dev/null +++ b/fmpcb_poly/test/t-revert_series_lagrange_fast.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("revert_series_lagrange_fast...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long qbits1, rbits1, rbits2, n; + fmpq_poly_t A, B; + fmpcb_poly_t a, b, c; + + qbits1 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + + do { + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_set_coeff_ui(A, 0, 0); + } while (A->length < 2 || fmpz_is_zero(A->coeffs + 1)); + + fmpq_poly_revert_series(B, A, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_revert_series_lagrange_fast(b, a, n, rbits2); + + if (!fmpcb_poly_contains_fmpq_poly(b, B)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits2 = %ld\n", n, rbits2); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(c, a); + fmpcb_poly_revert_series_lagrange_fast(c, c, n, rbits2); + if (!fmpcb_poly_equal(c, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/fmpcb_poly/test/t-revert_series_newton.c b/fmpcb_poly/test/t-revert_series_newton.c new file mode 100644 index 00000000..ae6cbd81 --- /dev/null +++ b/fmpcb_poly/test/t-revert_series_newton.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB 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. + + ARB 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 ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("revert_series_newton...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long qbits1, rbits1, rbits2, n; + fmpq_poly_t A, B; + fmpcb_poly_t a, b, c; + + qbits1 = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + n = 2 + n_randint(state, 25); + + fmpq_poly_init(A); + fmpq_poly_init(B); + + fmpcb_poly_init(a); + fmpcb_poly_init(b); + fmpcb_poly_init(c); + + do { + fmpq_poly_randtest(A, state, 1 + n_randint(state, 25), qbits1); + fmpq_poly_set_coeff_ui(A, 0, 0); + } while (A->length < 2 || fmpz_is_zero(A->coeffs + 1)); + + fmpq_poly_revert_series(B, A, n); + + fmpcb_poly_set_fmpq_poly(a, A, rbits1); + fmpcb_poly_revert_series_newton(b, a, n, rbits2); + + if (!fmpcb_poly_contains_fmpq_poly(b, B)) + { + printf("FAIL\n\n"); + printf("n = %ld, bits2 = %ld\n", n, rbits2); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("B = "); fmpq_poly_print(B); printf("\n\n"); + + printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n"); + + abort(); + } + + fmpcb_poly_set(c, a); + fmpcb_poly_revert_series_newton(c, c, n, rbits2); + if (!fmpcb_poly_equal(c, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + + fmpcb_poly_clear(a); + fmpcb_poly_clear(b); + fmpcb_poly_clear(c); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}