hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 17:31:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

implement Borel transform, binomial transform for fmprb_poly

Browse source
This commit is contained in:
Fredrik Johansson 2013-08-20 00:21:59 +02:00
parent f9c95b142a
commit 95f3462378
11 changed files with 879 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

51
doc/source/fmprb_poly.rst
View file

@ -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
-------------------------------------------------------------------------------

22
fmprb_poly.h
View file

@ -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,

63
fmprb_poly/binomial_transform.c Normal file
View file

@ -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);
}

85
fmprb_poly/binomial_transform_basecase.c Normal file
View file

@ -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);
}

96
fmprb_poly/binomial_transform_convolution.c Normal file
View file

@ -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 <math.h>
#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);
}

57
fmprb_poly/borel_transform.c Normal file
View file

@ -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);
}

57
fmprb_poly/inv_borel_transform.c Normal file
View file

@ -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);
}

119
fmprb_poly/test/t-binomial_transform.c Normal file
View file

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

119
fmprb_poly/test/t-binomial_transform_basecase.c Normal file
View file

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

119
fmprb_poly/test/t-binomial_transform_convolution.c Normal file
View file

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

91
fmprb_poly/test/t-borel_transform.c Normal file
View file

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