From 95f346237859ada3e4642789807235368b66de07 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 20 Aug 2013 00:21:59 +0200 Subject: [PATCH] implement Borel transform, binomial transform for fmprb_poly --- doc/source/fmprb_poly.rst | 51 ++++++++ fmprb_poly.h | 22 ++++ fmprb_poly/binomial_transform.c | 63 ++++++++++ fmprb_poly/binomial_transform_basecase.c | 85 +++++++++++++ fmprb_poly/binomial_transform_convolution.c | 96 ++++++++++++++ fmprb_poly/borel_transform.c | 57 +++++++++ fmprb_poly/inv_borel_transform.c | 57 +++++++++ fmprb_poly/test/t-binomial_transform.c | 119 ++++++++++++++++++ .../test/t-binomial_transform_basecase.c | 119 ++++++++++++++++++ .../test/t-binomial_transform_convolution.c | 119 ++++++++++++++++++ fmprb_poly/test/t-borel_transform.c | 91 ++++++++++++++ 11 files changed, 879 insertions(+) create mode 100644 fmprb_poly/binomial_transform.c create mode 100644 fmprb_poly/binomial_transform_basecase.c create mode 100644 fmprb_poly/binomial_transform_convolution.c create mode 100644 fmprb_poly/borel_transform.c create mode 100644 fmprb_poly/inv_borel_transform.c create mode 100644 fmprb_poly/test/t-binomial_transform.c create mode 100644 fmprb_poly/test/t-binomial_transform_basecase.c create mode 100644 fmprb_poly/test/t-binomial_transform_convolution.c create mode 100644 fmprb_poly/test/t-borel_transform.c diff --git a/doc/source/fmprb_poly.rst b/doc/source/fmprb_poly.rst index 5ff3827b..d5c44d89 100644 --- a/doc/source/fmprb_poly.rst +++ b/doc/source/fmprb_poly.rst @@ -527,6 +527,57 @@ Differentiation Sets *res* to the integral of *poly*. +Transforms +------------------------------------------------------------------------------- + +.. function:: void _fmprb_poly_borel_transform(fmprb_ptr res, fmprb_srcptr poly, long len, long prec) + +.. function:: void fmprb_poly_borel_transform(fmprb_poly_t res, const fmprb_poly_t poly, long prec) + + Computes the Borel transform of the input polynomial, mapping `\sum_k a_k x^k` + to `\sum_k (a_k / k!) x^k`. The underscore method allows aliasing. + +.. function:: void _fmprb_poly_inv_borel_transform(fmprb_ptr res, fmprb_srcptr poly, long len, long prec) + +.. function:: void fmprb_poly_inv_borel_transform(fmprb_poly_t res, const fmprb_poly_t poly, long prec) + + Computes the inverse Borel transform of the input polynomial, mapping `\sum_k a_k x^k` + to `\sum_k a_k k! x^k`. The underscore method allows aliasing. + +.. function:: void _fmprb_poly_binomial_transform_basecase(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec) + +.. function:: void fmprb_poly_binomial_transform_basecase(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec) + +.. function:: void _fmprb_poly_binomial_transform_convolution(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec) + +.. function:: void fmprb_poly_binomial_transform_convolution(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec) + +.. function:: void _fmprb_poly_binomial_transform(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec) + +.. function:: void fmprb_poly_binomial_transform(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec) + + Computes the binomial transform of the input truncated to length *len*. + The binomial transform maps the coefficients `a_k` in the input polynomial + to the coefficients `b_k` in the output polynomial via + `b_n = \sum_{k=0}^n (-1)^k {n \choose k} a_k`. + The binomial transform is equivalent to the power series composition + `f(x) \to (1-x)^{-1} f(x/(x-1))`, and is its own inverse. + + The *basecase* version evaluates coefficients one by one from the + definition, generating the binomial coefficients by a recurrence + relation. + + The *convolution* version uses the identity + `T(f(x)) = B^{-1}(e^x B(f(-x)))` where `T` denotes the binomial + transform operator and `B` denotes the Borel transform operator. + This only costs a single polynomial multiplication, plus some + scalar operations. + + The default version automatically chooses an algorithm. + + The underscore methods do not support aliasing, and assume that + the lengths are nonzero. + Special functions ------------------------------------------------------------------------------- diff --git a/fmprb_poly.h b/fmprb_poly.h index 0b3b5219..d2ae61d6 100644 --- a/fmprb_poly.h +++ b/fmprb_poly.h @@ -427,6 +427,28 @@ void _fmprb_poly_integral(fmprb_ptr res, fmprb_srcptr poly, long len, long prec) void fmprb_poly_integral(fmprb_poly_t res, const fmprb_poly_t poly, long prec); +/* Transforms */ + +void fmprb_poly_borel_transform(fmprb_poly_t res, const fmprb_poly_t poly, long prec); + +void _fmprb_poly_borel_transform(fmprb_ptr res, fmprb_srcptr poly, long len, long prec); + +void fmprb_poly_inv_borel_transform(fmprb_poly_t res, const fmprb_poly_t poly, long prec); + +void _fmprb_poly_inv_borel_transform(fmprb_ptr res, fmprb_srcptr poly, long len, long prec); + +void _fmprb_poly_binomial_transform_basecase(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec); + +void fmprb_poly_binomial_transform_basecase(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec); + +void _fmprb_poly_binomial_transform_convolution(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec); + +void fmprb_poly_binomial_transform_convolution(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec); + +void _fmprb_poly_binomial_transform(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec); + +void fmprb_poly_binomial_transform(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec); + /* Special functions */ void _fmprb_poly_rsqrt_series(fmprb_ptr g, diff --git a/fmprb_poly/binomial_transform.c b/fmprb_poly/binomial_transform.c new file mode 100644 index 00000000..ef167626 --- /dev/null +++ b/fmprb_poly/binomial_transform.c @@ -0,0 +1,63 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +void +_fmprb_poly_binomial_transform(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec) +{ + if (alen < 10 || len < 10) + _fmprb_poly_binomial_transform_basecase(b, a, alen, len, prec); + else + _fmprb_poly_binomial_transform_convolution(b, a, alen, len, prec); +} + +void +fmprb_poly_binomial_transform(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec) +{ + if (len == 0 || a->length == 0) + { + fmprb_poly_zero(b); + return; + } + + if (b == a) + { + fmprb_poly_t c; + fmprb_poly_init2(c, len); + _fmprb_poly_binomial_transform(c->coeffs, a->coeffs, a->length, len, prec); + fmprb_poly_swap(b, c); + fmprb_poly_clear(c); + } + else + { + fmprb_poly_fit_length(b, len); + _fmprb_poly_binomial_transform(b->coeffs, a->coeffs, a->length, len, prec); + } + + _fmprb_poly_set_length(b, len); + _fmprb_poly_normalise(b); +} + diff --git a/fmprb_poly/binomial_transform_basecase.c b/fmprb_poly/binomial_transform_basecase.c new file mode 100644 index 00000000..0d328a5d --- /dev/null +++ b/fmprb_poly/binomial_transform_basecase.c @@ -0,0 +1,85 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +void +_fmprb_poly_binomial_transform_basecase(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec) +{ + long n, k; + + fmpz_t t; + fmpz_init(t); + + for (n = 0; n < len; n++) + { + fmprb_zero(b + n); + + for (k = 0; k < FLINT_MIN(n + 1, alen); k++) + { + if (k == 0) + { + fmpz_one(t); + } + else + { + fmpz_mul_si(t, t, -(n - k + 1)); + fmpz_divexact_ui(t, t, k); + } + + fmprb_addmul_fmpz(b + n, a + k, t, prec); + } + } + + fmpz_clear(t); +} + +void +fmprb_poly_binomial_transform_basecase(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec) +{ + if (len == 0 || a->length == 0) + { + fmprb_poly_zero(b); + return; + } + + if (b == a) + { + fmprb_poly_t c; + fmprb_poly_init2(c, len); + _fmprb_poly_binomial_transform_basecase(c->coeffs, a->coeffs, a->length, len, prec); + fmprb_poly_swap(b, c); + fmprb_poly_clear(c); + } + else + { + fmprb_poly_fit_length(b, len); + _fmprb_poly_binomial_transform_basecase(b->coeffs, a->coeffs, a->length, len, prec); + } + + _fmprb_poly_set_length(b, len); + _fmprb_poly_normalise(b); +} + diff --git a/fmprb_poly/binomial_transform_convolution.c b/fmprb_poly/binomial_transform_convolution.c new file mode 100644 index 00000000..b1095392 --- /dev/null +++ b/fmprb_poly/binomial_transform_convolution.c @@ -0,0 +1,96 @@ +/*============================================================================= + + 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 +#include "fmprb_poly.h" + +void +_fmprb_poly_binomial_transform_convolution(fmprb_ptr b, fmprb_srcptr a, long alen, long len, long prec) +{ + long i, s; + fmprb_ptr c, d; + + alen = FLINT_MIN(alen, len); + + c = _fmprb_vec_init(alen); + d = _fmprb_vec_init(len); + + /* Hack: rescale x -> 2^s x such that (2^s)^n / n! ~= 1, + improving efficiency of block multiplication. The proper + solution is to add automatic near-optimal scaling in the + block multiplication. */ + s = (log(len) - 1) * 1.44269504088896 + 0.5; + + _fmprb_poly_borel_transform(c, a, alen, prec); + for (i = 1; i < alen; i += 2) + fmprb_neg(c + i, c + i); + + fmprb_one(d); + for (i = 1; i < len; i++) + fmprb_div_ui(d + i, d + i - 1, i, prec); + + for (i = 0; i < alen; i++) + fmprb_mul_2exp_si(c + i, c + i, i * s); + for (i = 0; i < len; i++) + fmprb_mul_2exp_si(d + i, d + i, i * s); + + _fmprb_poly_mullow(b, d, len, c, alen, len, prec); + + for (i = 0; i < len; i++) + fmprb_mul_2exp_si(b + i, b + i, -i * s); + + _fmprb_poly_inv_borel_transform(b, b, len, prec); + + _fmprb_vec_clear(c, alen); + _fmprb_vec_clear(d, len); +} + +void +fmprb_poly_binomial_transform_convolution(fmprb_poly_t b, const fmprb_poly_t a, long len, long prec) +{ + if (len == 0 || a->length == 0) + { + fmprb_poly_zero(b); + return; + } + + if (b == a) + { + fmprb_poly_t c; + fmprb_poly_init2(c, len); + _fmprb_poly_binomial_transform_convolution(c->coeffs, a->coeffs, a->length, len, prec); + fmprb_poly_swap(b, c); + fmprb_poly_clear(c); + } + else + { + fmprb_poly_fit_length(b, len); + _fmprb_poly_binomial_transform_convolution(b->coeffs, a->coeffs, a->length, len, prec); + } + + _fmprb_poly_set_length(b, len); + _fmprb_poly_normalise(b); +} + diff --git a/fmprb_poly/borel_transform.c b/fmprb_poly/borel_transform.c new file mode 100644 index 00000000..9ccc918d --- /dev/null +++ b/fmprb_poly/borel_transform.c @@ -0,0 +1,57 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +void +_fmprb_poly_borel_transform(fmprb_ptr res, fmprb_srcptr poly, long len, long prec) +{ + long i; + + fmprb_t t; + fmprb_init(t); + + fmprb_one(t); + + for (i = 0; i < len; i++) + { + if (i > 1) + fmprb_mul_ui(t, t, i, prec); + + fmprb_div(res + i, poly + i, t, prec); + } + + fmprb_clear(t); +} + +void +fmprb_poly_borel_transform(fmprb_poly_t res, const fmprb_poly_t poly, long prec) +{ + fmprb_poly_fit_length(res, poly->length); + _fmprb_poly_borel_transform(res->coeffs, poly->coeffs, poly->length, prec); + _fmprb_poly_set_length(res, poly->length); + _fmprb_poly_normalise(res); +} + diff --git a/fmprb_poly/inv_borel_transform.c b/fmprb_poly/inv_borel_transform.c new file mode 100644 index 00000000..28d0d55d --- /dev/null +++ b/fmprb_poly/inv_borel_transform.c @@ -0,0 +1,57 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +void +_fmprb_poly_inv_borel_transform(fmprb_ptr res, fmprb_srcptr poly, long len, long prec) +{ + long i; + + fmprb_t t; + fmprb_init(t); + + fmprb_one(t); + + for (i = 0; i < len; i++) + { + if (i > 1) + fmprb_mul_ui(t, t, i, prec); + + fmprb_mul(res + i, poly + i, t, prec); + } + + fmprb_clear(t); +} + +void +fmprb_poly_inv_borel_transform(fmprb_poly_t res, const fmprb_poly_t poly, long prec) +{ + fmprb_poly_fit_length(res, poly->length); + _fmprb_poly_inv_borel_transform(res->coeffs, poly->coeffs, poly->length, prec); + _fmprb_poly_set_length(res, poly->length); + _fmprb_poly_normalise(res); +} + diff --git a/fmprb_poly/test/t-binomial_transform.c b/fmprb_poly/test/t-binomial_transform.c new file mode 100644 index 00000000..7a921f7f --- /dev/null +++ b/fmprb_poly/test/t-binomial_transform.c @@ -0,0 +1,119 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("binomial_transform...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 5000; iter++) + { + fmprb_poly_t a, b, c, d; + long j, n, prec; + + fmprb_poly_init(a); + fmprb_poly_init(b); + fmprb_poly_init(c); + fmprb_poly_init(d); + + n = n_randint(state, 20); + prec = 2 + n_randint(state, 200); + + fmprb_poly_randtest(a, state, n, prec, 10); + fmprb_poly_randtest(b, state, n, prec, 10); + fmprb_poly_randtest(c, state, n, prec, 10); + + /* check self-inversion property */ + fmprb_poly_binomial_transform(b, a, n, prec); + fmprb_poly_binomial_transform(c, b, n, prec); + + fmprb_poly_set(d, a); + fmprb_poly_truncate(d, n); + + if (!fmprb_poly_contains(c, d)) + { + printf("FAIL (containment)\n\n"); + printf("n = %ld, prec = %ld\n\n", n, prec); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("c: "); fmprb_poly_printd(c, 15); printf("\n\n"); + printf("d: "); fmprb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + fmprb_poly_set(d, a); + fmprb_poly_binomial_transform(d, d, n, prec); + if (!fmprb_poly_equal(d, b)) + { + printf("FAIL (aliasing)\n\n"); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("d: "); fmprb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + /* compare with power series operations */ + fmprb_poly_zero(d); + for (j = 1; j < n; j++) + fmprb_poly_set_coeff_si(d, j, -1); + fmprb_poly_compose_series(c, a, d, n, prec); + for (j = 0; j < n; j++) + fmprb_poly_set_coeff_si(d, j, 1); + fmprb_poly_mullow(c, c, d, n, prec); + + if (!fmprb_poly_overlaps(b, c)) + { + printf("FAIL (power series)\n\n"); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("c: "); fmprb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmprb_poly_clear(a); + fmprb_poly_clear(b); + fmprb_poly_clear(c); + fmprb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/fmprb_poly/test/t-binomial_transform_basecase.c b/fmprb_poly/test/t-binomial_transform_basecase.c new file mode 100644 index 00000000..4a0d5a51 --- /dev/null +++ b/fmprb_poly/test/t-binomial_transform_basecase.c @@ -0,0 +1,119 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("binomial_transform_basecase...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 5000; iter++) + { + fmprb_poly_t a, b, c, d; + long j, n, prec; + + fmprb_poly_init(a); + fmprb_poly_init(b); + fmprb_poly_init(c); + fmprb_poly_init(d); + + n = n_randint(state, 20); + prec = 2 + n_randint(state, 200); + + fmprb_poly_randtest(a, state, n, prec, 10); + fmprb_poly_randtest(b, state, n, prec, 10); + fmprb_poly_randtest(c, state, n, prec, 10); + + /* check self-inversion property */ + fmprb_poly_binomial_transform_basecase(b, a, n, prec); + fmprb_poly_binomial_transform_basecase(c, b, n, prec); + + fmprb_poly_set(d, a); + fmprb_poly_truncate(d, n); + + if (!fmprb_poly_contains(c, d)) + { + printf("FAIL (containment)\n\n"); + printf("n = %ld, prec = %ld\n\n", n, prec); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("c: "); fmprb_poly_printd(c, 15); printf("\n\n"); + printf("d: "); fmprb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + fmprb_poly_set(d, a); + fmprb_poly_binomial_transform_basecase(d, d, n, prec); + if (!fmprb_poly_equal(d, b)) + { + printf("FAIL (aliasing)\n\n"); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("d: "); fmprb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + /* compare with power series operations */ + fmprb_poly_zero(d); + for (j = 1; j < n; j++) + fmprb_poly_set_coeff_si(d, j, -1); + fmprb_poly_compose_series(c, a, d, n, prec); + for (j = 0; j < n; j++) + fmprb_poly_set_coeff_si(d, j, 1); + fmprb_poly_mullow(c, c, d, n, prec); + + if (!fmprb_poly_overlaps(b, c)) + { + printf("FAIL (power series)\n\n"); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("c: "); fmprb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmprb_poly_clear(a); + fmprb_poly_clear(b); + fmprb_poly_clear(c); + fmprb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/fmprb_poly/test/t-binomial_transform_convolution.c b/fmprb_poly/test/t-binomial_transform_convolution.c new file mode 100644 index 00000000..69c63e2b --- /dev/null +++ b/fmprb_poly/test/t-binomial_transform_convolution.c @@ -0,0 +1,119 @@ +/*============================================================================= + + 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 "fmprb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("binomial_transform_convolution...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 5000; iter++) + { + fmprb_poly_t a, b, c, d; + long j, n, prec; + + fmprb_poly_init(a); + fmprb_poly_init(b); + fmprb_poly_init(c); + fmprb_poly_init(d); + + n = n_randint(state, 20); + prec = 2 + n_randint(state, 200); + + fmprb_poly_randtest(a, state, n, prec, 10); + fmprb_poly_randtest(b, state, n, prec, 10); + fmprb_poly_randtest(c, state, n, prec, 10); + + /* check self-inversion property */ + fmprb_poly_binomial_transform_convolution(b, a, n, prec); + fmprb_poly_binomial_transform_convolution(c, b, n, prec); + + fmprb_poly_set(d, a); + fmprb_poly_truncate(d, n); + + if (!fmprb_poly_contains(c, d)) + { + printf("FAIL (containment)\n\n"); + printf("n = %ld, prec = %ld\n\n", n, prec); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("c: "); fmprb_poly_printd(c, 15); printf("\n\n"); + printf("d: "); fmprb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + fmprb_poly_set(d, a); + fmprb_poly_binomial_transform_convolution(d, d, n, prec); + if (!fmprb_poly_equal(d, b)) + { + printf("FAIL (aliasing)\n\n"); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("d: "); fmprb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + /* compare with power series operations */ + fmprb_poly_zero(d); + for (j = 1; j < n; j++) + fmprb_poly_set_coeff_si(d, j, -1); + fmprb_poly_compose_series(c, a, d, n, prec); + for (j = 0; j < n; j++) + fmprb_poly_set_coeff_si(d, j, 1); + fmprb_poly_mullow(c, c, d, n, prec); + + if (!fmprb_poly_overlaps(b, c)) + { + printf("FAIL (power series)\n\n"); + + printf("a: "); fmprb_poly_printd(a, 15); printf("\n\n"); + printf("b: "); fmprb_poly_printd(b, 15); printf("\n\n"); + printf("c: "); fmprb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + fmprb_poly_clear(a); + fmprb_poly_clear(b); + fmprb_poly_clear(c); + fmprb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/fmprb_poly/test/t-borel_transform.c b/fmprb_poly/test/t-borel_transform.c new file mode 100644 index 00000000..23a2662e --- /dev/null +++ b/fmprb_poly/test/t-borel_transform.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 "fmprb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("borel_transform...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + fmprb_poly_t a, b, c, d; + long n, prec; + + fmprb_poly_init(a); + fmprb_poly_init(b); + fmprb_poly_init(c); + fmprb_poly_init(d); + + n = n_randint(state, 30); + prec = n_randint(state, 200); + + fmprb_poly_randtest(a, state, n, prec, 10); + fmprb_poly_randtest(b, state, n, prec, 10); + fmprb_poly_randtest(c, state, n, prec, 10); + + fmprb_poly_borel_transform(b, a, prec); + fmprb_poly_inv_borel_transform(c, b, prec); + + if (!fmprb_poly_contains(c, a)) + { + printf("FAIL (containment)\n\n"); + abort(); + } + + fmprb_poly_set(d, a); + fmprb_poly_borel_transform(d, d, prec); + if (!fmprb_poly_equal(d, b)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + fmprb_poly_set(d, b); + fmprb_poly_inv_borel_transform(d, d, prec); + if (!fmprb_poly_equal(d, c)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmprb_poly_clear(a); + fmprb_poly_clear(b); + fmprb_poly_clear(c); + fmprb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +