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

add fmpcb_poly compose, compose_series, revert methods

Browse source
This commit is contained in:
Fredrik Johansson 2013-07-27 11:37:26 +02:00
parent b45def0f92
commit e3b238884b
22 changed files with 2479 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

68
doc/source/fmpcb_poly.rst
View file

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

60
fmpcb_poly.h
View file

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

78
fmpcb_poly/compose.c Normal file
View file

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

196
fmpcb_poly/compose_divconquer.c Normal file
View file

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

125
fmpcb_poly/compose_horner.c Normal file
View file

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

112
fmpcb_poly/compose_series.c Normal file
View file

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

135
fmpcb_poly/compose_series_brent_kung.c Normal file
View file

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

120
fmpcb_poly/compose_series_horner.c Normal file
View file

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

91
fmpcb_poly/revert_series.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 "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);
}

122
fmpcb_poly/revert_series_lagrange.c Normal file
View file

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

147
fmpcb_poly/revert_series_lagrange_fast.c Normal file
View file

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

131
fmpcb_poly/revert_series_newton.c Normal file
View file

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

114
fmpcb_poly/test/t-compose.c Normal file
View file

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

114
fmpcb_poly/test/t-compose_divconquer.c Normal file
View file

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

114
fmpcb_poly/test/t-compose_horner.c Normal file
View file

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

116
fmpcb_poly/test/t-compose_series.c Normal file
View file

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

116
fmpcb_poly/test/t-compose_series_brent_kung.c Normal file
View file

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

116
fmpcb_poly/test/t-compose_series_horner.c Normal file
View file

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

101
fmpcb_poly/test/t-revert_series.c Normal file
View file

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

101
fmpcb_poly/test/t-revert_series_lagrange.c Normal file
View file

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

101
fmpcb_poly/test/t-revert_series_lagrange_fast.c Normal file
View file

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

101
fmpcb_poly/test/t-revert_series_newton.c Normal file
View file

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