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

don't require zero-padded input in series reversion

Browse source
This commit is contained in:
Fredrik Johansson 2014-01-15 11:40:11 +01:00
parent e6ecc103cf
commit 27ebe64238
12 changed files with 70 additions and 270 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
doc/source/fmpcb_poly.rst
View file

@ -350,19 +350,19 @@ Composition
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_ptr h, fmpcb_srcptr f, long flen, 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_ptr h, fmpcb_srcptr f, long flen, 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_ptr h, fmpcb_srcptr f, long flen, 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_ptr h, fmpcb_srcptr f, long flen, long n, long prec)
.. function:: void fmpcb_poly_revert_series(fmpcb_poly_t h, const fmpcb_poly_t f, long n, long prec)
@ -373,8 +373,8 @@ Composition
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.
linear term is nonzero. The underscore methods assume that *flen*
is at least 2, and do not support aliasing.
Evaluation
-------------------------------------------------------------------------------

12
doc/source/fmprb_poly.rst
View file

@ -365,19 +365,19 @@ Composition
with either input polynomial.
.. function:: void _fmprb_poly_revert_series_lagrange(fmprb_ptr h, fmprb_srcptr f, long n, long prec)
.. function:: void _fmprb_poly_revert_series_lagrange(fmprb_ptr h, fmprb_srcptr f, long flen, long n, long prec)
.. function:: void fmprb_poly_revert_series_lagrange(fmprb_poly_t h, const fmprb_poly_t f, long n, long prec)
.. function:: void _fmprb_poly_revert_series_newton(fmprb_ptr h, fmprb_srcptr f, long n, long prec)
.. function:: void _fmprb_poly_revert_series_newton(fmprb_ptr h, fmprb_srcptr f, long flen, long n, long prec)
.. function:: void fmprb_poly_revert_series_newton(fmprb_poly_t h, const fmprb_poly_t f, long n, long prec)
.. function:: void _fmprb_poly_revert_series_lagrange_fast(fmprb_ptr h, fmprb_srcptr f, long n, long prec)
.. function:: void _fmprb_poly_revert_series_lagrange_fast(fmprb_ptr h, fmprb_srcptr f, long flen, long n, long prec)
.. function:: void fmprb_poly_revert_series_lagrange_fast(fmprb_poly_t h, const fmprb_poly_t f, long n, long prec)
.. function:: void _fmprb_poly_revert_series(fmprb_ptr h, fmprb_srcptr f, long n, long prec)
.. function:: void _fmprb_poly_revert_series(fmprb_ptr h, fmprb_srcptr f, long flen, long n, long prec)
.. function:: void fmprb_poly_revert_series(fmprb_poly_t h, const fmprb_poly_t f, long n, long prec)
@ -388,8 +388,8 @@ Composition
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.
linear term is nonzero. The underscore methods assume that *flen*
is at least 2, and do not support aliasing.
Evaluation
-------------------------------------------------------------------------------

8
fmpcb_poly.h
View file

@ -335,16 +335,16 @@ void fmpcb_poly_compose_series(fmpcb_poly_t res,
/* Reversion */
void _fmpcb_poly_revert_series_lagrange(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec);
void _fmpcb_poly_revert_series_lagrange(fmpcb_ptr Qinv, fmpcb_srcptr Q, long Qlen, 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_ptr Qinv, fmpcb_srcptr Q, long Qlen, 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_ptr Qinv, fmpcb_srcptr Q, long Qlen, 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_ptr Qinv, fmpcb_srcptr Q, long Qlen, long n, long prec);
void fmpcb_poly_revert_series(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec);

37
fmpcb_poly/revert_series.c
View file

@ -27,65 +27,40 @@
void
_fmpcb_poly_revert_series(fmpcb_ptr Qinv,
fmpcb_srcptr Q, long n, long prec)
fmpcb_srcptr Q, long Qlen, long n, long prec)
{
_fmpcb_poly_revert_series_lagrange_fast(Qinv, Q, n, prec);
_fmpcb_poly_revert_series_lagrange_fast(Qinv, Q, Qlen, 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))
if (Qlen < 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);
_fmpcb_poly_revert_series(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmpcb_poly_t t;
fmpcb_poly_init2(t, n);
_fmpcb_poly_revert_series(t->coeffs, Qcopy, n, prec);
_fmpcb_poly_revert_series(t->coeffs, Q->coeffs, Qlen, 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);
}

37
fmpcb_poly/revert_series_lagrange.c
View file

@ -27,7 +27,7 @@
void
_fmpcb_poly_revert_series_lagrange(fmpcb_ptr Qinv,
fmpcb_srcptr Q, long n, long prec)
fmpcb_srcptr Q, long Qlen, long n, long prec)
{
long i;
fmpcb_ptr R, S, T, tmp;
@ -48,7 +48,7 @@ _fmpcb_poly_revert_series_lagrange(fmpcb_ptr Qinv,
fmpcb_zero(Qinv);
fmpcb_inv(Qinv + 1, Q + 1, prec);
_fmpcb_poly_inv_series(R, Q + 1, n - 1, n - 1, prec);
_fmpcb_poly_inv_series(R, Q + 1, FLINT_MIN(Qlen, n) - 1, n - 1, prec);
_fmpcb_vec_set(S, R, n - 1);
for (i = 2; i < n; i++)
@ -67,56 +67,31 @@ 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))
if (Qlen < 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);
_fmpcb_poly_revert_series_lagrange(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmpcb_poly_t t;
fmpcb_poly_init2(t, n);
_fmpcb_poly_revert_series_lagrange(t->coeffs, Qcopy, n, prec);
_fmpcb_poly_revert_series_lagrange(t->coeffs, Q->coeffs, Qlen, 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);
}

37
fmpcb_poly/revert_series_lagrange_fast.c
View file

@ -29,7 +29,7 @@
#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)
_fmpcb_poly_revert_series_lagrange_fast(fmpcb_ptr Qinv, fmpcb_srcptr Q, long Qlen, long n, long prec)
{
long i, j, k, m;
fmpcb_ptr R, S, T, tmp;
@ -54,7 +54,7 @@ _fmpcb_poly_revert_series_lagrange_fast(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n,
fmpcb_zero(Qinv);
fmpcb_inv(Qinv + 1, Q + 1, prec);
_fmpcb_poly_inv_series(Ri(1), Q + 1, n - 1, n - 1, prec);
_fmpcb_poly_inv_series(Ri(1), Q + 1, FLINT_MIN(Qlen, 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);
@ -92,56 +92,31 @@ 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))
if (Qlen < 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);
_fmpcb_poly_revert_series_lagrange_fast(Qinv->coeffs, Q->coeffs, Qlen, 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_revert_series_lagrange_fast(t->coeffs, Q->coeffs, Qlen, 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);
}

39
fmpcb_poly/revert_series_newton.c
View file

@ -28,7 +28,7 @@
#define CUTOFF 5
void
_fmpcb_poly_revert_series_newton(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long prec)
_fmpcb_poly_revert_series_newton(fmpcb_ptr Qinv, fmpcb_srcptr Q, long Qlen, long n, long prec)
{
long i, k, a[FLINT_BITS];
fmpcb_ptr T, U, V;
@ -52,13 +52,13 @@ _fmpcb_poly_revert_series_newton(fmpcb_ptr Qinv, fmpcb_srcptr Q, long n, long pr
while (k >= CUTOFF)
a[++i] = (k = (k + 1) / 2);
_fmpcb_poly_revert_series_lagrange(Qinv, Q, k, prec);
_fmpcb_poly_revert_series_lagrange(Qinv, Q, Qlen, 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_compose_series(T, Q, FLINT_MIN(Qlen, 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);
@ -76,56 +76,31 @@ 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))
if (Qlen < 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);
_fmpcb_poly_revert_series_newton(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmpcb_poly_t t;
fmpcb_poly_init2(t, n);
_fmpcb_poly_revert_series_newton(t->coeffs, Qcopy, n, prec);
_fmpcb_poly_revert_series_newton(t->coeffs, Q->coeffs, Qlen, 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);
}

8
fmprb_poly.h
View file

@ -333,16 +333,16 @@ void fmprb_poly_compose_series(fmprb_poly_t res,
/* Reversion */
void _fmprb_poly_revert_series_lagrange(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long prec);
void _fmprb_poly_revert_series_lagrange(fmprb_ptr Qinv, fmprb_srcptr Q, long Qlen, long n, long prec);
void fmprb_poly_revert_series_lagrange(fmprb_poly_t Qinv, const fmprb_poly_t Q, long n, long prec);
void _fmprb_poly_revert_series_newton(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long prec);
void _fmprb_poly_revert_series_newton(fmprb_ptr Qinv, fmprb_srcptr Q, long Qlen, long n, long prec);
void fmprb_poly_revert_series_newton(fmprb_poly_t Qinv, const fmprb_poly_t Q, long n, long prec);
void _fmprb_poly_revert_series_lagrange_fast(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long prec);
void _fmprb_poly_revert_series_lagrange_fast(fmprb_ptr Qinv, fmprb_srcptr Q, long Qlen, long n, long prec);
void fmprb_poly_revert_series_lagrange_fast(fmprb_poly_t Qinv, const fmprb_poly_t Q, long n, long prec);
void _fmprb_poly_revert_series(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long prec);
void _fmprb_poly_revert_series(fmprb_ptr Qinv, fmprb_srcptr Q, long Qlen, long n, long prec);
void fmprb_poly_revert_series(fmprb_poly_t Qinv, const fmprb_poly_t Q, long n, long prec);
/* Evaluation and interpolation */

37
fmprb_poly/revert_series.c
View file

@ -27,65 +27,40 @@
void
_fmprb_poly_revert_series(fmprb_ptr Qinv,
fmprb_srcptr Q, long n, long prec)
fmprb_srcptr Q, long Qlen, long n, long prec)
{
_fmprb_poly_revert_series_lagrange_fast(Qinv, Q, n, prec);
_fmprb_poly_revert_series_lagrange_fast(Qinv, Q, Qlen, n, prec);
}
void
fmprb_poly_revert_series(fmprb_poly_t Qinv,
const fmprb_poly_t Q, long n, long prec)
{
fmprb_ptr Qcopy;
int Qalloc;
long Qlen = Q->length;
if (Q->length < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
if (Qlen < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
{
printf("Exception (fmprb_poly_revert_series). Input must \n"
"have zero constant term and nonzero coefficient of x^1.\n");
abort();
}
if (n < 2)
{
fmprb_poly_zero(Qinv);
return;
}
if (Qlen >= n)
{
Qcopy = Q->coeffs;
Qalloc = 0;
}
else
{
long i;
Qcopy = _fmprb_vec_init(n);
for (i = 0; i < Qlen; i++)
Qcopy[i] = Q->coeffs[i];
Qalloc = 1;
}
if (Qinv != Q)
{
fmprb_poly_fit_length(Qinv, n);
_fmprb_poly_revert_series(Qinv->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmprb_poly_t t;
fmprb_poly_init2(t, n);
_fmprb_poly_revert_series(t->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series(t->coeffs, Q->coeffs, Qlen, n, prec);
fmprb_poly_swap(Qinv, t);
fmprb_poly_clear(t);
}
_fmprb_poly_set_length(Qinv, n);
_fmprb_poly_normalise(Qinv);
if (Qalloc)
flint_free(Qcopy);
}

37
fmprb_poly/revert_series_lagrange.c
View file

@ -27,7 +27,7 @@
void
_fmprb_poly_revert_series_lagrange(fmprb_ptr Qinv,
fmprb_srcptr Q, long n, long prec)
fmprb_srcptr Q, long Qlen, long n, long prec)
{
long i;
fmprb_ptr R, S, T, tmp;
@ -48,7 +48,7 @@ _fmprb_poly_revert_series_lagrange(fmprb_ptr Qinv,
fmprb_zero(Qinv);
fmprb_inv(Qinv + 1, Q + 1, prec);
_fmprb_poly_inv_series(R, Q + 1, n - 1, n - 1, prec);
_fmprb_poly_inv_series(R, Q + 1, FLINT_MIN(Qlen, n) - 1, n - 1, prec);
_fmprb_vec_set(S, R, n - 1);
for (i = 2; i < n; i++)
@ -67,56 +67,31 @@ void
fmprb_poly_revert_series_lagrange(fmprb_poly_t Qinv,
const fmprb_poly_t Q, long n, long prec)
{
fmprb_ptr Qcopy;
int Qalloc;
long Qlen = Q->length;
if (Q->length < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
if (Qlen < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
{
printf("Exception (fmprb_poly_revert_series_lagrange). Input must \n"
"have zero constant term and nonzero coefficient of x^1.\n");
abort();
}
if (n < 2)
{
fmprb_poly_zero(Qinv);
return;
}
if (Qlen >= n)
{
Qcopy = Q->coeffs;
Qalloc = 0;
}
else
{
long i;
Qcopy = _fmprb_vec_init(n);
for (i = 0; i < Qlen; i++)
Qcopy[i] = Q->coeffs[i];
Qalloc = 1;
}
if (Qinv != Q)
{
fmprb_poly_fit_length(Qinv, n);
_fmprb_poly_revert_series_lagrange(Qinv->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series_lagrange(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmprb_poly_t t;
fmprb_poly_init2(t, n);
_fmprb_poly_revert_series_lagrange(t->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series_lagrange(t->coeffs, Q->coeffs, Qlen, n, prec);
fmprb_poly_swap(Qinv, t);
fmprb_poly_clear(t);
}
_fmprb_poly_set_length(Qinv, n);
_fmprb_poly_normalise(Qinv);
if (Qalloc)
flint_free(Qcopy);
}

37
fmprb_poly/revert_series_lagrange_fast.c
View file

@ -29,7 +29,7 @@
#define Ri(ii) (R + (n-1)*((ii)-1))
void
_fmprb_poly_revert_series_lagrange_fast(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long prec)
_fmprb_poly_revert_series_lagrange_fast(fmprb_ptr Qinv, fmprb_srcptr Q, long Qlen, long n, long prec)
{
long i, j, k, m;
fmprb_ptr R, S, T, tmp;
@ -54,7 +54,7 @@ _fmprb_poly_revert_series_lagrange_fast(fmprb_ptr Qinv, fmprb_srcptr Q, long n,
fmprb_zero(Qinv);
fmprb_inv(Qinv + 1, Q + 1, prec);
_fmprb_poly_inv_series(Ri(1), Q + 1, n - 1, n - 1, prec);
_fmprb_poly_inv_series(Ri(1), Q + 1, FLINT_MIN(Qlen, n) - 1, n - 1, prec);
for (i = 2; i <= m; i++)
_fmprb_poly_mullow(Ri(i), Ri((i + 1) / 2), n - 1, Ri(i / 2), n - 1, n - 1, prec);
@ -92,56 +92,31 @@ void
fmprb_poly_revert_series_lagrange_fast(fmprb_poly_t Qinv,
const fmprb_poly_t Q, long n, long prec)
{
fmprb_ptr Qcopy;
int Qalloc;
long Qlen = Q->length;
if (Q->length < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
if (Qlen < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
{
printf("Exception (fmprb_poly_revert_series_lagrange_fast). Input \n"
"must have zero constant term and nonzero coefficient of x^1.\n");
abort();
}
if (n < 2)
{
fmprb_poly_zero(Qinv);
return;
}
if (Qlen >= n)
{
Qcopy = Q->coeffs;
Qalloc = 0;
}
else
{
long i;
Qcopy = _fmprb_vec_init(n);
for (i = 0; i < Qlen; i++)
Qcopy[i] = Q->coeffs[i];
Qalloc = 1;
}
if (Qinv != Q)
{
fmprb_poly_fit_length(Qinv, n);
_fmprb_poly_revert_series_lagrange_fast(Qinv->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series_lagrange_fast(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmprb_poly_t t;
fmprb_poly_init2(t, n);
_fmprb_poly_revert_series_lagrange_fast(t->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series_lagrange_fast(t->coeffs, Q->coeffs, Qlen, n, prec);
fmprb_poly_swap(Qinv, t);
fmprb_poly_clear(t);
}
_fmprb_poly_set_length(Qinv, n);
_fmprb_poly_normalise(Qinv);
if (Qalloc)
flint_free(Qcopy);
}

39
fmprb_poly/revert_series_newton.c
View file

@ -28,7 +28,7 @@
#define CUTOFF 5
void
_fmprb_poly_revert_series_newton(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long prec)
_fmprb_poly_revert_series_newton(fmprb_ptr Qinv, fmprb_srcptr Q, long Qlen, long n, long prec)
{
long i, k, a[FLINT_BITS];
fmprb_ptr T, U, V;
@ -52,13 +52,13 @@ _fmprb_poly_revert_series_newton(fmprb_ptr Qinv, fmprb_srcptr Q, long n, long pr
while (k >= CUTOFF)
a[++i] = (k = (k + 1) / 2);
_fmprb_poly_revert_series_lagrange(Qinv, Q, k, prec);
_fmprb_poly_revert_series_lagrange(Qinv, Q, Qlen, k, prec);
_fmprb_vec_zero(Qinv + k, n - k);
for (i--; i >= 0; i--)
{
k = a[i];
_fmprb_poly_compose_series(T, Q, k, Qinv, k, k, prec);
_fmprb_poly_compose_series(T, Q, FLINT_MIN(Qlen, k), Qinv, k, k, prec);
_fmprb_poly_derivative(U, T, k, prec); fmprb_zero(U + k - 1);
fmprb_zero(T + 1);
_fmprb_poly_div_series(V, T, k, U, k, k, prec);
@ -76,56 +76,31 @@ void
fmprb_poly_revert_series_newton(fmprb_poly_t Qinv,
const fmprb_poly_t Q, long n, long prec)
{
fmprb_ptr Qcopy;
int Qalloc;
long Qlen = Q->length;
if (Q->length < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
if (Qlen < 2 || !fmprb_is_zero(Q->coeffs)
|| fmprb_contains_zero(Q->coeffs + 1))
{
printf("Exception (fmprb_poly_revert_series_newton). Input must \n"
"have zero constant term and nonzero coefficient of x^1.\n");
abort();
}
if (n < 2)
{
fmprb_poly_zero(Qinv);
return;
}
if (Qlen >= n)
{
Qcopy = Q->coeffs;
Qalloc = 0;
}
else
{
long i;
Qcopy = _fmprb_vec_init(n);
for (i = 0; i < Qlen; i++)
Qcopy[i] = Q->coeffs[i];
Qalloc = 1;
}
if (Qinv != Q)
{
fmprb_poly_fit_length(Qinv, n);
_fmprb_poly_revert_series_newton(Qinv->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series_newton(Qinv->coeffs, Q->coeffs, Qlen, n, prec);
}
else
{
fmprb_poly_t t;
fmprb_poly_init2(t, n);
_fmprb_poly_revert_series_newton(t->coeffs, Qcopy, n, prec);
_fmprb_poly_revert_series_newton(t->coeffs, Q->coeffs, Qlen, n, prec);
fmprb_poly_swap(Qinv, t);
fmprb_poly_clear(t);
}
_fmprb_poly_set_length(Qinv, n);
_fmprb_poly_normalise(Qinv);
if (Qalloc)
flint_free(Qcopy);
}