From ed2990d5f8b4f79aba783ac3390fb0cb5777c6b9 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 10 Feb 2015 17:45:20 +0100 Subject: [PATCH] add sin_pi_series and friends --- acb_poly.h | 22 ++++- acb_poly/cos_pi_series.c | 82 ++++++++++++++++ acb_poly/gamma_series.c | 6 +- acb_poly/rgamma_series.c | 6 +- acb_poly/sin_cos_pi_series.c | 90 +++++++++++++++++ acb_poly/sin_cos_series.c | 4 +- acb_poly/sin_cos_series_basecase.c | 19 +++- acb_poly/sin_cos_series_tangent.c | 22 ++++- acb_poly/sin_pi_series.c | 75 ++++++++++++++ acb_poly/tan_series.c | 2 +- acb_poly/test/t-cos_pi_series.c | 101 +++++++++++++++++++ acb_poly/test/t-sin_cos_pi_series.c | 114 ++++++++++++++++++++++ acb_poly/test/t-sin_cos_series_basecase.c | 8 +- acb_poly/test/t-sin_cos_series_tangent.c | 8 +- acb_poly/test/t-sin_pi_series.c | 101 +++++++++++++++++++ acb_poly/zeta_series.c | 10 +- arb_poly.h | 22 ++++- arb_poly/cos_pi_series.c | 82 ++++++++++++++++ arb_poly/gamma_series.c | 6 +- arb_poly/rgamma_series.c | 6 +- arb_poly/sin_cos_pi_series.c | 90 +++++++++++++++++ arb_poly/sin_cos_series.c | 4 +- arb_poly/sin_cos_series_basecase.c | 19 +++- arb_poly/sin_cos_series_tangent.c | 22 ++++- arb_poly/sin_pi_series.c | 75 ++++++++++++++ arb_poly/tan_series.c | 2 +- arb_poly/test/t-acos_series.c | 2 +- arb_poly/test/t-asin_series.c | 2 +- arb_poly/test/t-cos_pi_series.c | 101 +++++++++++++++++++ arb_poly/test/t-sin_cos_pi_series.c | 114 ++++++++++++++++++++++ arb_poly/test/t-sin_cos_series_basecase.c | 8 +- arb_poly/test/t-sin_cos_series_tangent.c | 8 +- arb_poly/test/t-sin_pi_series.c | 101 +++++++++++++++++++ arb_poly/zeta_series.c | 10 +- doc/source/acb_poly.rst | 26 ++++- doc/source/arb_poly.rst | 26 ++++- 36 files changed, 1318 insertions(+), 78 deletions(-) create mode 100644 acb_poly/cos_pi_series.c create mode 100644 acb_poly/sin_cos_pi_series.c create mode 100644 acb_poly/sin_pi_series.c create mode 100644 acb_poly/test/t-cos_pi_series.c create mode 100644 acb_poly/test/t-sin_cos_pi_series.c create mode 100644 acb_poly/test/t-sin_pi_series.c create mode 100644 arb_poly/cos_pi_series.c create mode 100644 arb_poly/sin_cos_pi_series.c create mode 100644 arb_poly/sin_pi_series.c create mode 100644 arb_poly/test/t-cos_pi_series.c create mode 100644 arb_poly/test/t-sin_cos_pi_series.c create mode 100644 arb_poly/test/t-sin_pi_series.c diff --git a/acb_poly.h b/acb_poly.h index febed996..193a480a 100644 --- a/acb_poly.h +++ b/acb_poly.h @@ -493,16 +493,16 @@ void _acb_poly_exp_series(acb_ptr f, acb_srcptr h, long hlen, long n, long prec) void acb_poly_exp_series(acb_poly_t f, const acb_poly_t h, long n, long prec); void _acb_poly_sin_cos_series_basecase(acb_ptr s, - acb_ptr c, acb_srcptr h, long hlen, long n, long prec); + acb_ptr c, acb_srcptr h, long hlen, long n, long prec, int times_pi); void acb_poly_sin_cos_series_basecase(acb_poly_t s, acb_poly_t c, - const acb_poly_t h, long n, long prec); + const acb_poly_t h, long n, long prec, int times_pi); void _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, - const acb_srcptr h, long hlen, long len, long prec); + const acb_srcptr h, long hlen, long len, long prec, int times_pi); void acb_poly_sin_cos_series_tangent(acb_poly_t s, acb_poly_t c, - const acb_poly_t h, long n, long prec); + const acb_poly_t h, long n, long prec, int times_pi); void _acb_poly_sin_cos_series(acb_ptr s, acb_ptr c, const acb_srcptr h, long hlen, long len, long prec); @@ -518,6 +518,20 @@ void _acb_poly_cos_series(acb_ptr g, acb_srcptr h, long hlen, long n, long prec) void acb_poly_cos_series(acb_poly_t g, const acb_poly_t h, long n, long prec); +void _acb_poly_sin_cos_pi_series(acb_ptr s, acb_ptr c, + const acb_srcptr h, long hlen, long len, long prec); + +void acb_poly_sin_cos_pi_series(acb_poly_t s, acb_poly_t c, + const acb_poly_t h, long n, long prec); + +void _acb_poly_sin_pi_series(acb_ptr g, acb_srcptr h, long hlen, long n, long prec); + +void acb_poly_sin_pi_series(acb_poly_t g, const acb_poly_t h, long n, long prec); + +void _acb_poly_cos_pi_series(acb_ptr g, acb_srcptr h, long hlen, long n, long prec); + +void acb_poly_cos_pi_series(acb_poly_t g, const acb_poly_t h, long n, long prec); + void _acb_poly_tan_series(acb_ptr g, acb_srcptr h, long hlen, long len, long prec); void acb_poly_tan_series(acb_poly_t g, const acb_poly_t h, long n, long prec); diff --git a/acb_poly/cos_pi_series.c b/acb_poly/cos_pi_series.c new file mode 100644 index 00000000..2cfde5c0 --- /dev/null +++ b/acb_poly/cos_pi_series.c @@ -0,0 +1,82 @@ +/*============================================================================= + + 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 "acb_poly.h" + +void +_acb_poly_cos_pi_series(acb_ptr g, acb_srcptr h, long hlen, long n, long prec) +{ + hlen = FLINT_MIN(hlen, n); + + if (hlen == 1) + { + acb_cos_pi(g, h, prec); + _acb_vec_zero(g + 1, n - 1); + } + else if (n == 2) + { + acb_t t; + acb_init(t); + acb_sin_cos_pi(t, g, h, prec); + acb_neg(t, t); + acb_mul(g + 1, h + 1, t, prec); /* safe since hlen >= 2 */ + acb_const_pi(t, prec); + acb_mul(g + 1, g + 1, t, prec); + acb_clear(t); + } + else + { + acb_ptr t = _acb_vec_init(n); + _acb_poly_sin_cos_pi_series(t, g, h, hlen, n, prec); + _acb_vec_clear(t, n); + } +} + +void +acb_poly_cos_pi_series(acb_poly_t g, const acb_poly_t h, long n, long prec) +{ + long hlen = h->length; + + if (n == 0) + { + acb_poly_zero(g); + return; + } + + if (hlen == 0) + { + acb_poly_one(g); + return; + } + + if (hlen == 1) + n = 1; + + acb_poly_fit_length(g, n); + _acb_poly_cos_pi_series(g->coeffs, h->coeffs, hlen, n, prec); + _acb_poly_set_length(g, n); + _acb_poly_normalise(g); +} + diff --git a/acb_poly/gamma_series.c b/acb_poly/gamma_series.c index afed8b83..05ac42f7 100644 --- a/acb_poly/gamma_series.c +++ b/acb_poly/gamma_series.c @@ -221,9 +221,9 @@ _acb_poly_gamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long prec acb_neg(u + i, u + i); /* v = 1/sin(pi x) */ - acb_const_pi(f + 1, wp); - acb_mul(f, h, f + 1, wp); - _acb_poly_sin_series(t, f, 2, len, wp); + acb_set(f, h); + acb_one(f + 1); + _acb_poly_sin_pi_series(t, f, 2, len, wp); _acb_poly_inv_series(v, t, len, len, wp); _acb_poly_mullow(t, u, len, v, len, len, wp); diff --git a/acb_poly/rgamma_series.c b/acb_poly/rgamma_series.c index 9ac03d17..8dc6aad0 100644 --- a/acb_poly/rgamma_series.c +++ b/acb_poly/rgamma_series.c @@ -84,9 +84,9 @@ _acb_poly_rgamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long pre acb_neg(u + i, u + i); /* v = sin(pi x) */ - acb_const_pi(f + 1, wp); - acb_mul(f, h, f + 1, wp); - _acb_poly_sin_series(v, f, 2, len, wp); + acb_set(f, h); + acb_one(f + 1); + _acb_poly_sin_pi_series(v, f, 2, len, wp); _acb_poly_mullow(t, u, len, v, len, len, wp); diff --git a/acb_poly/sin_cos_pi_series.c b/acb_poly/sin_cos_pi_series.c new file mode 100644 index 00000000..186e25b8 --- /dev/null +++ b/acb_poly/sin_cos_pi_series.c @@ -0,0 +1,90 @@ +/*============================================================================= + + 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 "acb_poly.h" + +#define TANGENT_CUTOFF 80 + +void +_acb_poly_sin_cos_pi_series(acb_ptr s, acb_ptr c, const acb_srcptr h, long hlen, long n, long prec) +{ + hlen = FLINT_MIN(hlen, n); + + if (hlen == 1) + { + acb_sin_cos_pi(s, c, h, prec); + _acb_vec_zero(s + 1, n - 1); + _acb_vec_zero(c + 1, n - 1); + } + else if (n == 2) + { + acb_t t; + acb_init(t); + acb_const_pi(t, prec); + acb_mul(t, t, h + 1, prec); + acb_sin_cos_pi(s, c, h, prec); + acb_mul(s + 1, c, t, prec); + acb_neg(t, t); + acb_mul(c + 1, s, t, prec); + acb_clear(t); + } + else if (hlen < TANGENT_CUTOFF) + _acb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 1); + else + _acb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 1); +} + +void +acb_poly_sin_cos_pi_series(acb_poly_t s, acb_poly_t c, + const acb_poly_t h, long n, long prec) +{ + long hlen = h->length; + + if (n == 0) + { + acb_poly_zero(s); + acb_poly_zero(c); + return; + } + + if (hlen == 0) + { + acb_poly_zero(s); + acb_poly_one(c); + return; + } + + if (hlen == 1) + n = 1; + + acb_poly_fit_length(s, n); + acb_poly_fit_length(c, n); + _acb_poly_sin_cos_pi_series(s->coeffs, c->coeffs, h->coeffs, hlen, n, prec); + _acb_poly_set_length(s, n); + _acb_poly_normalise(s); + _acb_poly_set_length(c, n); + _acb_poly_normalise(c); +} + diff --git a/acb_poly/sin_cos_series.c b/acb_poly/sin_cos_series.c index d5002a71..5721034f 100644 --- a/acb_poly/sin_cos_series.c +++ b/acb_poly/sin_cos_series.c @@ -50,9 +50,9 @@ _acb_poly_sin_cos_series(acb_ptr s, acb_ptr c, const acb_srcptr h, long hlen, lo acb_clear(t); } else if (hlen < TANGENT_CUTOFF) - _acb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec); + _acb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 0); else - _acb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec); + _acb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 0); } void diff --git a/acb_poly/sin_cos_series_basecase.c b/acb_poly/sin_cos_series_basecase.c index 33f1d250..59fbb085 100644 --- a/acb_poly/sin_cos_series_basecase.c +++ b/acb_poly/sin_cos_series_basecase.c @@ -26,14 +26,17 @@ #include "acb_poly.h" void -_acb_poly_sin_cos_series_basecase(acb_ptr s, - acb_ptr c, acb_srcptr h, long hlen, long n, long prec) +_acb_poly_sin_cos_series_basecase(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, + long n, long prec, int times_pi) { long j, k, alen = FLINT_MIN(n, hlen); acb_ptr a; acb_t t, u; - acb_sin_cos(s, c, h, prec); + if (times_pi) + acb_sin_cos_pi(s, c, h, prec); + else + acb_sin_cos(s, c, h, prec); if (hlen == 1) { @@ -49,6 +52,12 @@ _acb_poly_sin_cos_series_basecase(acb_ptr s, for (k = 1; k < alen; k++) acb_mul_ui(a + k, h + k, k, prec); + if (times_pi) + { + acb_const_pi(t, prec); + _acb_vec_scalar_mul(a + 1, a + 1, alen - 1, t, prec); + } + for (k = 1; k < n; k++) { acb_zero(t); @@ -71,7 +80,7 @@ _acb_poly_sin_cos_series_basecase(acb_ptr s, void acb_poly_sin_cos_series_basecase(acb_poly_t s, acb_poly_t c, - const acb_poly_t h, long n, long prec) + const acb_poly_t h, long n, long prec, int times_pi) { long hlen = h->length; @@ -91,7 +100,7 @@ acb_poly_sin_cos_series_basecase(acb_poly_t s, acb_poly_t c, acb_poly_fit_length(s, n); acb_poly_fit_length(c, n); - _acb_poly_sin_cos_series_basecase(s->coeffs, c->coeffs, h->coeffs, hlen, n, prec); + _acb_poly_sin_cos_series_basecase(s->coeffs, c->coeffs, h->coeffs, hlen, n, prec, times_pi); _acb_poly_set_length(s, n); _acb_poly_normalise(s); _acb_poly_set_length(c, n); diff --git a/acb_poly/sin_cos_series_tangent.c b/acb_poly/sin_cos_series_tangent.c index 37b3b34a..c8f22a11 100644 --- a/acb_poly/sin_cos_series_tangent.c +++ b/acb_poly/sin_cos_series_tangent.c @@ -27,7 +27,7 @@ void _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, - const acb_srcptr h, long hlen, long len, long prec) + const acb_srcptr h, long hlen, long len, long prec, int times_pi) { acb_ptr t, u, v; acb_t s0, c0; @@ -35,7 +35,10 @@ _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, if (hlen == 1) { - acb_sin_cos(s, c, h, prec); + if (times_pi) + acb_sin_cos_pi(s, c, h, prec); + else + acb_sin_cos(s, c, h, prec); _acb_vec_zero(s + 1, len - 1); _acb_vec_zero(c + 1, len - 1); return; @@ -54,11 +57,20 @@ _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, v = u + len; /* sin, cos of h0 */ - acb_sin_cos(s0, c0, h, prec); + if (times_pi) + acb_sin_cos_pi(s0, c0, h, prec); + else + acb_sin_cos(s0, c0, h, prec); /* t = tan((h-h0)/2) */ acb_zero(u); _acb_vec_scalar_mul_2exp_si(u + 1, h + 1, hlen - 1, -1); + if (times_pi) + { + acb_const_pi(t, prec); + _acb_vec_scalar_mul(u + 1, u + 1, hlen - 1, t, prec); + } + _acb_poly_tan_series(t, u, hlen, len, prec); /* v = 1 + t^2 */ @@ -97,7 +109,7 @@ _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, void acb_poly_sin_cos_series_tangent(acb_poly_t s, acb_poly_t c, - const acb_poly_t h, long n, long prec) + const acb_poly_t h, long n, long prec, int times_pi) { long hlen = h->length; @@ -118,7 +130,7 @@ acb_poly_sin_cos_series_tangent(acb_poly_t s, acb_poly_t c, acb_poly_fit_length(s, n); acb_poly_fit_length(c, n); _acb_poly_sin_cos_series_tangent(s->coeffs, c->coeffs, - h->coeffs, hlen, n, prec); + h->coeffs, hlen, n, prec, times_pi); _acb_poly_set_length(s, n); _acb_poly_normalise(s); _acb_poly_set_length(c, n); diff --git a/acb_poly/sin_pi_series.c b/acb_poly/sin_pi_series.c new file mode 100644 index 00000000..27d580da --- /dev/null +++ b/acb_poly/sin_pi_series.c @@ -0,0 +1,75 @@ +/*============================================================================= + + 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 "acb_poly.h" + +void +_acb_poly_sin_pi_series(acb_ptr g, acb_srcptr h, long hlen, long n, long prec) +{ + hlen = FLINT_MIN(hlen, n); + + if (hlen == 1) + { + acb_sin_pi(g, h, prec); + _acb_vec_zero(g + 1, n - 1); + } + else if (n == 2) + { + acb_t t; + acb_init(t); + acb_sin_cos_pi(g, t, h, prec); + acb_mul(g + 1, h + 1, t, prec); /* safe since hlen >= 2 */ + acb_const_pi(t, prec); + acb_mul(g + 1, g + 1, t, prec); + acb_clear(t); + } + else + { + acb_ptr t = _acb_vec_init(n); + _acb_poly_sin_cos_pi_series(g, t, h, hlen, n, prec); + _acb_vec_clear(t, n); + } +} + +void +acb_poly_sin_pi_series(acb_poly_t g, const acb_poly_t h, long n, long prec) +{ + long hlen = h->length; + + if (hlen == 0 || n == 0) + { + acb_poly_zero(g); + return; + } + + if (hlen == 1) + n = 1; + + acb_poly_fit_length(g, n); + _acb_poly_sin_pi_series(g->coeffs, h->coeffs, hlen, n, prec); + _acb_poly_set_length(g, n); + _acb_poly_normalise(g); +} + diff --git a/acb_poly/tan_series.c b/acb_poly/tan_series.c index 5ae443bf..1ba2a839 100644 --- a/acb_poly/tan_series.c +++ b/acb_poly/tan_series.c @@ -58,7 +58,7 @@ _acb_poly_tan_series(acb_ptr g, NEWTON_INIT(TAN_NEWTON_CUTOFF, len) NEWTON_BASECASE(n) - _acb_poly_sin_cos_series_basecase(t, u, h, hlen, n, prec); + _acb_poly_sin_cos_series_basecase(t, u, h, hlen, n, prec, 0); _acb_poly_div_series(g, t, n, u, n, n, prec); NEWTON_END_BASECASE diff --git a/acb_poly/test/t-cos_pi_series.c b/acb_poly/test/t-cos_pi_series.c new file mode 100644 index 00000000..266bbfa0 --- /dev/null +++ b/acb_poly/test/t-cos_pi_series.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("cos_pi_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, bits1, bits2, bits3; + acb_poly_t S, A, B, C; + acb_t pi; + + bits1 = 2 + n_randint(state, 200); + bits2 = 2 + n_randint(state, 200); + bits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + acb_poly_init(S); + acb_poly_init(A); + acb_poly_init(B); + acb_poly_init(C); + acb_init(pi); + + acb_poly_randtest(S, state, m, bits1, 3); + acb_poly_randtest(A, state, m, bits1, 3); + acb_poly_randtest(B, state, m, bits1, 3); + + acb_poly_cos_pi_series(A, S, n1, bits2); + + acb_const_pi(pi, bits3); + acb_poly_set_acb(B, pi); + acb_poly_mul(B, S, B, bits3); + acb_poly_cos_series(B, B, n2, bits3); + + acb_poly_set(C, A); + acb_poly_truncate(C, FLINT_MIN(n1, n2)); + acb_poly_truncate(B, FLINT_MIN(n1, n2)); + + if (!acb_poly_overlaps(B, C)) + { + printf("FAIL\n\n"); + printf("S = "); acb_poly_printd(S, 15); printf("\n\n"); + printf("A = "); acb_poly_printd(A, 15); printf("\n\n"); + printf("B = "); acb_poly_printd(B, 15); printf("\n\n"); + abort(); + } + + acb_poly_cos_pi_series(S, S, n1, bits2); + + if (!acb_poly_overlaps(A, S)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + acb_poly_clear(S); + acb_poly_clear(A); + acb_poly_clear(B); + acb_poly_clear(C); + acb_clear(pi); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_poly/test/t-sin_cos_pi_series.c b/acb_poly/test/t-sin_cos_pi_series.c new file mode 100644 index 00000000..5a3f7ba9 --- /dev/null +++ b/acb_poly/test/t-sin_cos_pi_series.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("sin_cos_pi_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n, qbits, rbits1, rbits2; + fmpq_poly_t A, B; + acb_poly_t a, b, c, d, e; + + qbits = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n = 1 + n_randint(state, 30); + + fmpq_poly_init(A); + fmpq_poly_init(B); + acb_poly_init(a); + acb_poly_init(b); + acb_poly_init(c); + acb_poly_init(d); + acb_poly_init(e); + + fmpq_poly_randtest(A, state, m, qbits); + acb_poly_set_fmpq_poly(a, A, rbits1); + + acb_poly_sin_cos_pi_series(b, c, a, n, rbits2); + + /* Check sin(x)^2 + cos(x)^2 = 1 */ + acb_poly_mullow(d, b, b, n, rbits2); + acb_poly_mullow(e, c, c, n, rbits2); + acb_poly_add(d, d, e, rbits2); + + fmpq_poly_one(B); + if (!acb_poly_contains_fmpq_poly(d, B)) + { + printf("FAIL\n\n"); + printf("bits2 = %ld\n", rbits2); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); acb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + acb_poly_set(d, a); + acb_poly_sin_cos_pi_series(d, c, d, n, rbits2); + if (!acb_poly_equal(b, d)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + acb_poly_set(d, a); + acb_poly_sin_cos_pi_series(b, d, d, n, rbits2); + if (!acb_poly_equal(c, d)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + acb_poly_clear(a); + acb_poly_clear(b); + acb_poly_clear(c); + acb_poly_clear(d); + acb_poly_clear(e); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb_poly/test/t-sin_cos_series_basecase.c b/acb_poly/test/t-sin_cos_series_basecase.c index 82e3b90a..08e235d2 100644 --- a/acb_poly/test/t-sin_cos_series_basecase.c +++ b/acb_poly/test/t-sin_cos_series_basecase.c @@ -40,12 +40,14 @@ int main() long m, n, rbits1, rbits2; fmpq_poly_t B; acb_poly_t a, b, c, d, e; + int times_pi; rbits1 = 2 + n_randint(state, 200); rbits2 = 2 + n_randint(state, 200); m = 1 + n_randint(state, 30); n = 1 + n_randint(state, 30); + times_pi = n_randint(state, 2); fmpq_poly_init(B); acb_poly_init(a); @@ -56,7 +58,7 @@ int main() acb_poly_randtest(a, state, m, rbits1, 10); - acb_poly_sin_cos_series_basecase(b, c, a, n, rbits2); + acb_poly_sin_cos_series_basecase(b, c, a, n, rbits2, times_pi); /* Check sin(x)^2 + cos(x)^2 = 1 */ acb_poly_mullow(d, b, b, n, rbits2); @@ -78,7 +80,7 @@ int main() } acb_poly_set(d, a); - acb_poly_sin_cos_series_basecase(d, c, d, n, rbits2); + acb_poly_sin_cos_series_basecase(d, c, d, n, rbits2, times_pi); if (!acb_poly_equal(b, d)) { printf("FAIL (aliasing 1)\n\n"); @@ -86,7 +88,7 @@ int main() } acb_poly_set(d, a); - acb_poly_sin_cos_series_basecase(b, d, d, n, rbits2); + acb_poly_sin_cos_series_basecase(b, d, d, n, rbits2, times_pi); if (!acb_poly_equal(c, d)) { printf("FAIL (aliasing 2)\n\n"); diff --git a/acb_poly/test/t-sin_cos_series_tangent.c b/acb_poly/test/t-sin_cos_series_tangent.c index 1ca602b0..782e8569 100644 --- a/acb_poly/test/t-sin_cos_series_tangent.c +++ b/acb_poly/test/t-sin_cos_series_tangent.c @@ -40,12 +40,14 @@ int main() long m, n, rbits1, rbits2; fmpq_poly_t B; acb_poly_t a, b, c, d, e; + int times_pi; rbits1 = 2 + n_randint(state, 200); rbits2 = 2 + n_randint(state, 200); m = 1 + n_randint(state, 30); n = 1 + n_randint(state, 30); + times_pi = n_randint(state, 2); fmpq_poly_init(B); acb_poly_init(a); @@ -56,7 +58,7 @@ int main() acb_poly_randtest(a, state, m, rbits1, 10); - acb_poly_sin_cos_series_tangent(b, c, a, n, rbits2); + acb_poly_sin_cos_series_tangent(b, c, a, n, rbits2, times_pi); /* Check sin(x)^2 + cos(x)^2 = 1 */ acb_poly_mullow(d, b, b, n, rbits2); @@ -78,7 +80,7 @@ int main() } acb_poly_set(d, a); - acb_poly_sin_cos_series_tangent(d, c, d, n, rbits2); + acb_poly_sin_cos_series_tangent(d, c, d, n, rbits2, times_pi); if (!acb_poly_equal(b, d)) { printf("FAIL (aliasing 1)\n\n"); @@ -86,7 +88,7 @@ int main() } acb_poly_set(d, a); - acb_poly_sin_cos_series_tangent(b, d, d, n, rbits2); + acb_poly_sin_cos_series_tangent(b, d, d, n, rbits2, times_pi); if (!acb_poly_equal(c, d)) { printf("FAIL (aliasing 2)\n\n"); diff --git a/acb_poly/test/t-sin_pi_series.c b/acb_poly/test/t-sin_pi_series.c new file mode 100644 index 00000000..c5bc543e --- /dev/null +++ b/acb_poly/test/t-sin_pi_series.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("sin_pi_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, bits1, bits2, bits3; + acb_poly_t S, A, B, C; + acb_t pi; + + bits1 = 2 + n_randint(state, 200); + bits2 = 2 + n_randint(state, 200); + bits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + acb_poly_init(S); + acb_poly_init(A); + acb_poly_init(B); + acb_poly_init(C); + acb_init(pi); + + acb_poly_randtest(S, state, m, bits1, 3); + acb_poly_randtest(A, state, m, bits1, 3); + acb_poly_randtest(B, state, m, bits1, 3); + + acb_poly_sin_pi_series(A, S, n1, bits2); + + acb_const_pi(pi, bits3); + acb_poly_set_acb(B, pi); + acb_poly_mul(B, S, B, bits3); + acb_poly_sin_series(B, B, n2, bits3); + + acb_poly_set(C, A); + acb_poly_truncate(C, FLINT_MIN(n1, n2)); + acb_poly_truncate(B, FLINT_MIN(n1, n2)); + + if (!acb_poly_overlaps(B, C)) + { + printf("FAIL\n\n"); + printf("S = "); acb_poly_printd(S, 15); printf("\n\n"); + printf("A = "); acb_poly_printd(A, 15); printf("\n\n"); + printf("B = "); acb_poly_printd(B, 15); printf("\n\n"); + abort(); + } + + acb_poly_sin_pi_series(S, S, n1, bits2); + + if (!acb_poly_overlaps(A, S)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + acb_poly_clear(S); + acb_poly_clear(A); + acb_poly_clear(B); + acb_poly_clear(C); + acb_clear(pi); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_poly/zeta_series.c b/acb_poly/zeta_series.c index a7bd7770..9c679d4a 100644 --- a/acb_poly/zeta_series.c +++ b/acb_poly/zeta_series.c @@ -94,11 +94,11 @@ _acb_poly_zeta_series(acb_ptr res, acb_srcptr h, long hlen, const acb_t a, int d acb_mul_2exp_si(pi, pi, -1); /* s2 = sin(pi*s/2) / pi */ - acb_mul_2exp_si(pi, pi, -1); - acb_mul(f, pi, h, prec); - acb_set(f + 1, pi); - acb_mul_2exp_si(pi, pi, 1); - _acb_poly_sin_series(s2, f, 2, len, prec); + acb_set(f, h); + acb_one(f + 1); + acb_mul_2exp_si(f, f, -1); + acb_mul_2exp_si(f + 1, f + 1, -1); + _acb_poly_sin_pi_series(s2, f, 2, len, prec); _acb_vec_scalar_div(s2, s2, len, pi, prec); /* s3 = gamma(1-s) */ diff --git a/arb_poly.h b/arb_poly.h index 7ba83041..42967b44 100644 --- a/arb_poly.h +++ b/arb_poly.h @@ -509,16 +509,16 @@ void _arb_poly_exp_series(arb_ptr f, arb_srcptr h, long hlen, long n, long prec) void arb_poly_exp_series(arb_poly_t f, const arb_poly_t h, long n, long prec); void _arb_poly_sin_cos_series_basecase(arb_ptr s, - arb_ptr c, arb_srcptr h, long hlen, long n, long prec); + arb_ptr c, arb_srcptr h, long hlen, long n, long prec, int times_pi); void arb_poly_sin_cos_series_basecase(arb_poly_t s, arb_poly_t c, - const arb_poly_t h, long n, long prec); + const arb_poly_t h, long n, long prec, int times_pi); void _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, - const arb_srcptr h, long hlen, long len, long prec); + const arb_srcptr h, long hlen, long len, long prec, int times_pi); void arb_poly_sin_cos_series_tangent(arb_poly_t s, arb_poly_t c, - const arb_poly_t h, long n, long prec); + const arb_poly_t h, long n, long prec, int times_pi); void _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, const arb_srcptr h, long hlen, long len, long prec); @@ -526,6 +526,12 @@ void _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, void arb_poly_sin_cos_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, long n, long prec); +void _arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, + const arb_srcptr h, long hlen, long len, long prec); + +void arb_poly_sin_cos_pi_series(arb_poly_t s, arb_poly_t c, + const arb_poly_t h, long n, long prec); + void _arb_poly_sin_series(arb_ptr g, arb_srcptr h, long hlen, long n, long prec); void arb_poly_sin_series(arb_poly_t g, const arb_poly_t h, long n, long prec); @@ -534,6 +540,14 @@ void _arb_poly_cos_series(arb_ptr g, arb_srcptr h, long hlen, long n, long prec) void arb_poly_cos_series(arb_poly_t g, const arb_poly_t h, long n, long prec); +void _arb_poly_sin_pi_series(arb_ptr g, arb_srcptr h, long hlen, long n, long prec); + +void arb_poly_sin_pi_series(arb_poly_t g, const arb_poly_t h, long n, long prec); + +void _arb_poly_cos_pi_series(arb_ptr g, arb_srcptr h, long hlen, long n, long prec); + +void arb_poly_cos_pi_series(arb_poly_t g, const arb_poly_t h, long n, long prec); + void _arb_poly_tan_series(arb_ptr g, arb_srcptr h, long hlen, long len, long prec); void arb_poly_tan_series(arb_poly_t g, const arb_poly_t h, long n, long prec); diff --git a/arb_poly/cos_pi_series.c b/arb_poly/cos_pi_series.c new file mode 100644 index 00000000..6ba7c894 --- /dev/null +++ b/arb_poly/cos_pi_series.c @@ -0,0 +1,82 @@ +/*============================================================================= + + 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 "arb_poly.h" + +void +_arb_poly_cos_pi_series(arb_ptr g, arb_srcptr h, long hlen, long n, long prec) +{ + hlen = FLINT_MIN(hlen, n); + + if (hlen == 1) + { + arb_cos_pi(g, h, prec); + _arb_vec_zero(g + 1, n - 1); + } + else if (n == 2) + { + arb_t t; + arb_init(t); + arb_sin_cos_pi(t, g, h, prec); + arb_neg(t, t); + arb_mul(g + 1, h + 1, t, prec); /* safe since hlen >= 2 */ + arb_const_pi(t, prec); + arb_mul(g + 1, g + 1, t, prec); + arb_clear(t); + } + else + { + arb_ptr t = _arb_vec_init(n); + _arb_poly_sin_cos_pi_series(t, g, h, hlen, n, prec); + _arb_vec_clear(t, n); + } +} + +void +arb_poly_cos_pi_series(arb_poly_t g, const arb_poly_t h, long n, long prec) +{ + long hlen = h->length; + + if (n == 0) + { + arb_poly_zero(g); + return; + } + + if (hlen == 0) + { + arb_poly_one(g); + return; + } + + if (hlen == 1) + n = 1; + + arb_poly_fit_length(g, n); + _arb_poly_cos_pi_series(g->coeffs, h->coeffs, hlen, n, prec); + _arb_poly_set_length(g, n); + _arb_poly_normalise(g); +} + diff --git a/arb_poly/gamma_series.c b/arb_poly/gamma_series.c index 15533193..3f4d22e0 100644 --- a/arb_poly/gamma_series.c +++ b/arb_poly/gamma_series.c @@ -266,9 +266,9 @@ _arb_poly_gamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec arb_neg(u + i, u + i); /* v = 1/sin(pi x) */ - arb_const_pi(f + 1, wp); - arb_mul(f, h, f + 1, wp); - _arb_poly_sin_series(t, f, 2, len, wp); + arb_set(f, h); + arb_one(f + 1); + _arb_poly_sin_pi_series(t, f, 2, len, wp); _arb_poly_inv_series(v, t, len, len, wp); _arb_poly_mullow(t, u, len, v, len, len, wp); diff --git a/arb_poly/rgamma_series.c b/arb_poly/rgamma_series.c index ed4f689c..8f65cac1 100644 --- a/arb_poly/rgamma_series.c +++ b/arb_poly/rgamma_series.c @@ -129,9 +129,9 @@ _arb_poly_rgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long pre arb_neg(u + i, u + i); /* v = sin(pi x) */ - arb_const_pi(f + 1, wp); - arb_mul(f, h, f + 1, wp); - _arb_poly_sin_series(v, f, 2, len, wp); + arb_set(f, h); + arb_one(f + 1); + _arb_poly_sin_pi_series(v, f, 2, len, wp); _arb_poly_mullow(t, u, len, v, len, len, wp); diff --git a/arb_poly/sin_cos_pi_series.c b/arb_poly/sin_cos_pi_series.c new file mode 100644 index 00000000..75b7de77 --- /dev/null +++ b/arb_poly/sin_cos_pi_series.c @@ -0,0 +1,90 @@ +/*============================================================================= + + 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 "arb_poly.h" + +#define TANGENT_CUTOFF 240 + +void +_arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, const arb_srcptr h, long hlen, long n, long prec) +{ + hlen = FLINT_MIN(hlen, n); + + if (hlen == 1) + { + arb_sin_cos_pi(s, c, h, prec); + _arb_vec_zero(s + 1, n - 1); + _arb_vec_zero(c + 1, n - 1); + } + else if (n == 2) + { + arb_t t; + arb_init(t); + arb_const_pi(t, prec); + arb_mul(t, t, h + 1, prec); + arb_sin_cos_pi(s, c, h, prec); + arb_mul(s + 1, c, t, prec); + arb_neg(t, t); + arb_mul(c + 1, s, t, prec); + arb_clear(t); + } + else if (hlen < TANGENT_CUTOFF) + _arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 1); + else + _arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 1); +} + +void +arb_poly_sin_cos_pi_series(arb_poly_t s, arb_poly_t c, + const arb_poly_t h, long n, long prec) +{ + long hlen = h->length; + + if (n == 0) + { + arb_poly_zero(s); + arb_poly_zero(c); + return; + } + + if (hlen == 0) + { + arb_poly_zero(s); + arb_poly_one(c); + return; + } + + if (hlen == 1) + n = 1; + + arb_poly_fit_length(s, n); + arb_poly_fit_length(c, n); + _arb_poly_sin_cos_pi_series(s->coeffs, c->coeffs, h->coeffs, hlen, n, prec); + _arb_poly_set_length(s, n); + _arb_poly_normalise(s); + _arb_poly_set_length(c, n); + _arb_poly_normalise(c); +} + diff --git a/arb_poly/sin_cos_series.c b/arb_poly/sin_cos_series.c index fa8abc30..ae169afa 100644 --- a/arb_poly/sin_cos_series.c +++ b/arb_poly/sin_cos_series.c @@ -50,9 +50,9 @@ _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, const arb_srcptr h, long hlen, lo arb_clear(t); } else if (hlen < TANGENT_CUTOFF) - _arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec); + _arb_poly_sin_cos_series_basecase(s, c, h, hlen, n, prec, 0); else - _arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec); + _arb_poly_sin_cos_series_tangent(s, c, h, hlen, n, prec, 0); } void diff --git a/arb_poly/sin_cos_series_basecase.c b/arb_poly/sin_cos_series_basecase.c index 68a513f5..a21a663a 100644 --- a/arb_poly/sin_cos_series_basecase.c +++ b/arb_poly/sin_cos_series_basecase.c @@ -26,14 +26,17 @@ #include "arb_poly.h" void -_arb_poly_sin_cos_series_basecase(arb_ptr s, - arb_ptr c, arb_srcptr h, long hlen, long n, long prec) +_arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, + long n, long prec, int times_pi) { long j, k, alen = FLINT_MIN(n, hlen); arb_ptr a; arb_t t, u; - arb_sin_cos(s, c, h, prec); + if (times_pi) + arb_sin_cos_pi(s, c, h, prec); + else + arb_sin_cos(s, c, h, prec); if (hlen == 1) { @@ -49,6 +52,12 @@ _arb_poly_sin_cos_series_basecase(arb_ptr s, for (k = 1; k < alen; k++) arb_mul_ui(a + k, h + k, k, prec); + if (times_pi) + { + arb_const_pi(t, prec); + _arb_vec_scalar_mul(a + 1, a + 1, alen - 1, t, prec); + } + for (k = 1; k < n; k++) { arb_zero(t); @@ -71,7 +80,7 @@ _arb_poly_sin_cos_series_basecase(arb_ptr s, void arb_poly_sin_cos_series_basecase(arb_poly_t s, arb_poly_t c, - const arb_poly_t h, long n, long prec) + const arb_poly_t h, long n, long prec, int times_pi) { long hlen = h->length; @@ -91,7 +100,7 @@ arb_poly_sin_cos_series_basecase(arb_poly_t s, arb_poly_t c, arb_poly_fit_length(s, n); arb_poly_fit_length(c, n); - _arb_poly_sin_cos_series_basecase(s->coeffs, c->coeffs, h->coeffs, hlen, n, prec); + _arb_poly_sin_cos_series_basecase(s->coeffs, c->coeffs, h->coeffs, hlen, n, prec, times_pi); _arb_poly_set_length(s, n); _arb_poly_normalise(s); _arb_poly_set_length(c, n); diff --git a/arb_poly/sin_cos_series_tangent.c b/arb_poly/sin_cos_series_tangent.c index c9e1f46f..999cafdf 100644 --- a/arb_poly/sin_cos_series_tangent.c +++ b/arb_poly/sin_cos_series_tangent.c @@ -27,7 +27,7 @@ void _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, - const arb_srcptr h, long hlen, long len, long prec) + const arb_srcptr h, long hlen, long len, long prec, int times_pi) { arb_ptr t, u, v; arb_t s0, c0; @@ -35,7 +35,10 @@ _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, if (hlen == 1) { - arb_sin_cos(s, c, h, prec); + if (times_pi) + arb_sin_cos_pi(s, c, h, prec); + else + arb_sin_cos(s, c, h, prec); _arb_vec_zero(s + 1, len - 1); _arb_vec_zero(c + 1, len - 1); return; @@ -54,11 +57,20 @@ _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, v = u + len; /* sin, cos of h0 */ - arb_sin_cos(s0, c0, h, prec); + if (times_pi) + arb_sin_cos_pi(s0, c0, h, prec); + else + arb_sin_cos(s0, c0, h, prec); /* t = tan((h-h0)/2) */ arb_zero(u); _arb_vec_scalar_mul_2exp_si(u + 1, h + 1, hlen - 1, -1); + if (times_pi) + { + arb_const_pi(t, prec); + _arb_vec_scalar_mul(u + 1, u + 1, hlen - 1, t, prec); + } + _arb_poly_tan_series(t, u, hlen, len, prec); /* v = 1 + t^2 */ @@ -97,7 +109,7 @@ _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, void arb_poly_sin_cos_series_tangent(arb_poly_t s, arb_poly_t c, - const arb_poly_t h, long n, long prec) + const arb_poly_t h, long n, long prec, int times_pi) { long hlen = h->length; @@ -118,7 +130,7 @@ arb_poly_sin_cos_series_tangent(arb_poly_t s, arb_poly_t c, arb_poly_fit_length(s, n); arb_poly_fit_length(c, n); _arb_poly_sin_cos_series_tangent(s->coeffs, c->coeffs, - h->coeffs, hlen, n, prec); + h->coeffs, hlen, n, prec, times_pi); _arb_poly_set_length(s, n); _arb_poly_normalise(s); _arb_poly_set_length(c, n); diff --git a/arb_poly/sin_pi_series.c b/arb_poly/sin_pi_series.c new file mode 100644 index 00000000..ff788cf6 --- /dev/null +++ b/arb_poly/sin_pi_series.c @@ -0,0 +1,75 @@ +/*============================================================================= + + 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 "arb_poly.h" + +void +_arb_poly_sin_pi_series(arb_ptr g, arb_srcptr h, long hlen, long n, long prec) +{ + hlen = FLINT_MIN(hlen, n); + + if (hlen == 1) + { + arb_sin_pi(g, h, prec); + _arb_vec_zero(g + 1, n - 1); + } + else if (n == 2) + { + arb_t t; + arb_init(t); + arb_sin_cos_pi(g, t, h, prec); + arb_mul(g + 1, h + 1, t, prec); /* safe since hlen >= 2 */ + arb_const_pi(t, prec); + arb_mul(g + 1, g + 1, t, prec); + arb_clear(t); + } + else + { + arb_ptr t = _arb_vec_init(n); + _arb_poly_sin_cos_pi_series(g, t, h, hlen, n, prec); + _arb_vec_clear(t, n); + } +} + +void +arb_poly_sin_pi_series(arb_poly_t g, const arb_poly_t h, long n, long prec) +{ + long hlen = h->length; + + if (hlen == 0 || n == 0) + { + arb_poly_zero(g); + return; + } + + if (hlen == 1) + n = 1; + + arb_poly_fit_length(g, n); + _arb_poly_sin_pi_series(g->coeffs, h->coeffs, hlen, n, prec); + _arb_poly_set_length(g, n); + _arb_poly_normalise(g); +} + diff --git a/arb_poly/tan_series.c b/arb_poly/tan_series.c index c0afd01f..6a0498e6 100644 --- a/arb_poly/tan_series.c +++ b/arb_poly/tan_series.c @@ -58,7 +58,7 @@ _arb_poly_tan_series(arb_ptr g, NEWTON_INIT(TAN_NEWTON_CUTOFF, len) NEWTON_BASECASE(n) - _arb_poly_sin_cos_series_basecase(t, u, h, hlen, n, prec); + _arb_poly_sin_cos_series_basecase(t, u, h, hlen, n, prec, 0); _arb_poly_div_series(g, t, n, u, n, n, prec); NEWTON_END_BASECASE diff --git a/arb_poly/test/t-acos_series.c b/arb_poly/test/t-acos_series.c index 2f634d6a..20d1afa4 100644 --- a/arb_poly/test/t-acos_series.c +++ b/arb_poly/test/t-acos_series.c @@ -62,7 +62,7 @@ int main() arb_poly_acos_series(b, a, n, rbits2); /* Check cos(acos(x)) = x */ - arb_poly_sin_cos_series_basecase(d, c, b, n, rbits2); + arb_poly_sin_cos_series_basecase(d, c, b, n, rbits2, 0); fmpq_poly_truncate(A, n); if (!arb_poly_contains_fmpq_poly(c, A)) diff --git a/arb_poly/test/t-asin_series.c b/arb_poly/test/t-asin_series.c index 9dab9ece..a28f6432 100644 --- a/arb_poly/test/t-asin_series.c +++ b/arb_poly/test/t-asin_series.c @@ -62,7 +62,7 @@ int main() arb_poly_asin_series(b, a, n, rbits2); /* Check sin(asin(x)) = x */ - arb_poly_sin_cos_series_basecase(c, d, b, n, rbits2); + arb_poly_sin_cos_series_basecase(c, d, b, n, rbits2, 0); fmpq_poly_truncate(A, n); if (!arb_poly_contains_fmpq_poly(c, A)) diff --git a/arb_poly/test/t-cos_pi_series.c b/arb_poly/test/t-cos_pi_series.c new file mode 100644 index 00000000..2b2c1ba2 --- /dev/null +++ b/arb_poly/test/t-cos_pi_series.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("cos_pi_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, bits1, bits2, bits3; + arb_poly_t S, A, B, C; + arb_t pi; + + bits1 = 2 + n_randint(state, 200); + bits2 = 2 + n_randint(state, 200); + bits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + arb_poly_init(S); + arb_poly_init(A); + arb_poly_init(B); + arb_poly_init(C); + arb_init(pi); + + arb_poly_randtest(S, state, m, bits1, 3); + arb_poly_randtest(A, state, m, bits1, 3); + arb_poly_randtest(B, state, m, bits1, 3); + + arb_poly_cos_pi_series(A, S, n1, bits2); + + arb_const_pi(pi, bits3); + arb_poly_set_arb(B, pi); + arb_poly_mul(B, S, B, bits3); + arb_poly_cos_series(B, B, n2, bits3); + + arb_poly_set(C, A); + arb_poly_truncate(C, FLINT_MIN(n1, n2)); + arb_poly_truncate(B, FLINT_MIN(n1, n2)); + + if (!arb_poly_overlaps(B, C)) + { + printf("FAIL\n\n"); + printf("S = "); arb_poly_printd(S, 15); printf("\n\n"); + printf("A = "); arb_poly_printd(A, 15); printf("\n\n"); + printf("B = "); arb_poly_printd(B, 15); printf("\n\n"); + abort(); + } + + arb_poly_cos_pi_series(S, S, n1, bits2); + + if (!arb_poly_overlaps(A, S)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + arb_poly_clear(S); + arb_poly_clear(A); + arb_poly_clear(B); + arb_poly_clear(C); + arb_clear(pi); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb_poly/test/t-sin_cos_pi_series.c b/arb_poly/test/t-sin_cos_pi_series.c new file mode 100644 index 00000000..4b2736d1 --- /dev/null +++ b/arb_poly/test/t-sin_cos_pi_series.c @@ -0,0 +1,114 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("sin_cos_pi_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n, qbits, rbits1, rbits2; + fmpq_poly_t A, B; + arb_poly_t a, b, c, d, e; + + qbits = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n = 1 + n_randint(state, 30); + + fmpq_poly_init(A); + fmpq_poly_init(B); + arb_poly_init(a); + arb_poly_init(b); + arb_poly_init(c); + arb_poly_init(d); + arb_poly_init(e); + + fmpq_poly_randtest(A, state, m, qbits); + arb_poly_set_fmpq_poly(a, A, rbits1); + + arb_poly_sin_cos_pi_series(b, c, a, n, rbits2); + + /* Check sin(x)^2 + cos(x)^2 = 1 */ + arb_poly_mullow(d, b, b, n, rbits2); + arb_poly_mullow(e, c, c, n, rbits2); + arb_poly_add(d, d, e, rbits2); + + fmpq_poly_one(B); + if (!arb_poly_contains_fmpq_poly(d, B)) + { + printf("FAIL\n\n"); + printf("bits2 = %ld\n", rbits2); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); arb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + arb_poly_set(d, a); + arb_poly_sin_cos_pi_series(d, c, d, n, rbits2); + if (!arb_poly_equal(b, d)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + + arb_poly_set(d, a); + arb_poly_sin_cos_pi_series(b, d, d, n, rbits2); + if (!arb_poly_equal(c, d)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + fmpq_poly_clear(B); + arb_poly_clear(a); + arb_poly_clear(b); + arb_poly_clear(c); + arb_poly_clear(d); + arb_poly_clear(e); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb_poly/test/t-sin_cos_series_basecase.c b/arb_poly/test/t-sin_cos_series_basecase.c index 9e671d7d..23d864ab 100644 --- a/arb_poly/test/t-sin_cos_series_basecase.c +++ b/arb_poly/test/t-sin_cos_series_basecase.c @@ -40,6 +40,7 @@ int main() long m, n, qbits, rbits1, rbits2; fmpq_poly_t A, B; arb_poly_t a, b, c, d, e; + int times_pi; qbits = 2 + n_randint(state, 200); rbits1 = 2 + n_randint(state, 200); @@ -47,6 +48,7 @@ int main() m = 1 + n_randint(state, 30); n = 1 + n_randint(state, 30); + times_pi = n_randint(state, 2); fmpq_poly_init(A); fmpq_poly_init(B); @@ -59,7 +61,7 @@ int main() fmpq_poly_randtest(A, state, m, qbits); arb_poly_set_fmpq_poly(a, A, rbits1); - arb_poly_sin_cos_series_basecase(b, c, a, n, rbits2); + arb_poly_sin_cos_series_basecase(b, c, a, n, rbits2, times_pi); /* Check sin(x)^2 + cos(x)^2 = 1 */ arb_poly_mullow(d, b, b, n, rbits2); @@ -82,7 +84,7 @@ int main() } arb_poly_set(d, a); - arb_poly_sin_cos_series_basecase(d, c, d, n, rbits2); + arb_poly_sin_cos_series_basecase(d, c, d, n, rbits2, times_pi); if (!arb_poly_equal(b, d)) { printf("FAIL (aliasing 1)\n\n"); @@ -90,7 +92,7 @@ int main() } arb_poly_set(d, a); - arb_poly_sin_cos_series_basecase(b, d, d, n, rbits2); + arb_poly_sin_cos_series_basecase(b, d, d, n, rbits2, times_pi); if (!arb_poly_equal(c, d)) { printf("FAIL (aliasing 2)\n\n"); diff --git a/arb_poly/test/t-sin_cos_series_tangent.c b/arb_poly/test/t-sin_cos_series_tangent.c index 8526dfab..a6214445 100644 --- a/arb_poly/test/t-sin_cos_series_tangent.c +++ b/arb_poly/test/t-sin_cos_series_tangent.c @@ -40,6 +40,7 @@ int main() long m, n, qbits, rbits1, rbits2; fmpq_poly_t A, B; arb_poly_t a, b, c, d, e; + int times_pi; qbits = 2 + n_randint(state, 200); rbits1 = 2 + n_randint(state, 200); @@ -47,6 +48,7 @@ int main() m = 1 + n_randint(state, 30); n = 1 + n_randint(state, 30); + times_pi = n_randint(state, 2); fmpq_poly_init(A); fmpq_poly_init(B); @@ -59,7 +61,7 @@ int main() fmpq_poly_randtest(A, state, m, qbits); arb_poly_set_fmpq_poly(a, A, rbits1); - arb_poly_sin_cos_series_tangent(b, c, a, n, rbits2); + arb_poly_sin_cos_series_tangent(b, c, a, n, rbits2, times_pi); /* Check sin(x)^2 + cos(x)^2 = 1 */ arb_poly_mullow(d, b, b, n, rbits2); @@ -82,7 +84,7 @@ int main() } arb_poly_set(d, a); - arb_poly_sin_cos_series_tangent(d, c, d, n, rbits2); + arb_poly_sin_cos_series_tangent(d, c, d, n, rbits2, times_pi); if (!arb_poly_equal(b, d)) { printf("FAIL (aliasing 1)\n\n"); @@ -90,7 +92,7 @@ int main() } arb_poly_set(d, a); - arb_poly_sin_cos_series_tangent(b, d, d, n, rbits2); + arb_poly_sin_cos_series_tangent(b, d, d, n, rbits2, times_pi); if (!arb_poly_equal(c, d)) { printf("FAIL (aliasing 2)\n\n"); diff --git a/arb_poly/test/t-sin_pi_series.c b/arb_poly/test/t-sin_pi_series.c new file mode 100644 index 00000000..998efbc1 --- /dev/null +++ b/arb_poly/test/t-sin_pi_series.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("sin_pi_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, bits1, bits2, bits3; + arb_poly_t S, A, B, C; + arb_t pi; + + bits1 = 2 + n_randint(state, 200); + bits2 = 2 + n_randint(state, 200); + bits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + arb_poly_init(S); + arb_poly_init(A); + arb_poly_init(B); + arb_poly_init(C); + arb_init(pi); + + arb_poly_randtest(S, state, m, bits1, 3); + arb_poly_randtest(A, state, m, bits1, 3); + arb_poly_randtest(B, state, m, bits1, 3); + + arb_poly_sin_pi_series(A, S, n1, bits2); + + arb_const_pi(pi, bits3); + arb_poly_set_arb(B, pi); + arb_poly_mul(B, S, B, bits3); + arb_poly_sin_series(B, B, n2, bits3); + + arb_poly_set(C, A); + arb_poly_truncate(C, FLINT_MIN(n1, n2)); + arb_poly_truncate(B, FLINT_MIN(n1, n2)); + + if (!arb_poly_overlaps(B, C)) + { + printf("FAIL\n\n"); + printf("S = "); arb_poly_printd(S, 15); printf("\n\n"); + printf("A = "); arb_poly_printd(A, 15); printf("\n\n"); + printf("B = "); arb_poly_printd(B, 15); printf("\n\n"); + abort(); + } + + arb_poly_sin_pi_series(S, S, n1, bits2); + + if (!arb_poly_overlaps(A, S)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + arb_poly_clear(S); + arb_poly_clear(A); + arb_poly_clear(B); + arb_poly_clear(C); + arb_clear(pi); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb_poly/zeta_series.c b/arb_poly/zeta_series.c index 6c7784c4..1d11c5b3 100644 --- a/arb_poly/zeta_series.c +++ b/arb_poly/zeta_series.c @@ -100,11 +100,11 @@ _arb_poly_zeta_series(arb_ptr res, arb_srcptr h, long hlen, const arb_t a, int d arb_mul_2exp_si(pi, pi, -1); /* s2 = sin(pi*s/2) / pi */ - arb_mul_2exp_si(pi, pi, -1); - arb_mul(f, pi, h, prec); - arb_set(f + 1, pi); - arb_mul_2exp_si(pi, pi, 1); - _arb_poly_sin_series(s2, f, 2, len, prec); + arb_set(f, h); + arb_one(f + 1); + arb_mul_2exp_si(f, f, -1); + arb_mul_2exp_si(f + 1, f + 1, -1); + _arb_poly_sin_pi_series(s2, f, 2, len, prec); _arb_vec_scalar_div(s2, s2, len, pi, prec); /* s3 = gamma(1-s) */ diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst index 0f064c6e..2f960fb1 100644 --- a/doc/source/acb_poly.rst +++ b/doc/source/acb_poly.rst @@ -670,13 +670,13 @@ Elementary functions The underscore methods support aliasing and allow the input to be shorter than the output, but require the lengths to be nonzero. -.. function:: void _acb_poly_sin_cos_series_basecase(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, long n, long prec) +.. function:: void _acb_poly_sin_cos_series_basecase(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, long n, long prec, int times_pi) -.. function:: void acb_poly_sin_cos_series_basecase(acb_poly_t s, acb_poly_t c, const acb_poly_t h, long n, long prec) +.. function:: void acb_poly_sin_cos_series_basecase(acb_poly_t s, acb_poly_t c, const acb_poly_t h, long n, long prec, int times_pi) -.. function:: void _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, long n, long prec) +.. function:: void _acb_poly_sin_cos_series_tangent(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, long n, long prec, int times_pi) -.. function:: void acb_poly_sin_cos_series_tangent(acb_poly_t s, acb_poly_t c, const acb_poly_t h, long n, long prec) +.. function:: void acb_poly_sin_cos_series_tangent(acb_poly_t s, acb_poly_t c, const acb_poly_t h, long n, long prec, int times_pi) .. function:: void _acb_poly_sin_cos_series(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, long n, long prec) @@ -701,6 +701,9 @@ Elementary functions The default version automatically selects between the *basecase* and *tangent* algorithms depending on the input. + The *basecase* and *tangent* versions take a flag *times_pi* + specifying that the input is to be multiplied by `\pi`. + The underscore methods support aliasing and require the lengths to be nonzero. .. function:: void _acb_poly_sin_series(acb_ptr s, acb_srcptr h, long hlen, long n, long prec) @@ -715,6 +718,21 @@ Elementary functions simply wrap :func:`_acb_poly_sin_cos_series`. The underscore methods support aliasing and require the lengths to be nonzero. +.. function:: void _acb_poly_sin_cos_pi_series(acb_ptr s, acb_ptr c, acb_srcptr h, long hlen, long n, long prec) + +.. function:: void acb_poly_sin_cos_pi_series(acb_poly_t s, acb_poly_t c, const acb_poly_t h, long n, long prec) + +.. function:: void _acb_poly_sin_pi_series(acb_ptr s, acb_srcptr h, long hlen, long n, long prec) + +.. function:: void acb_poly_sin_pi_series(acb_poly_t s, const acb_poly_t h, long n, long prec) + +.. function:: void _acb_poly_cos_pi_series(acb_ptr c, acb_srcptr h, long hlen, long n, long prec) + +.. function:: void acb_poly_cos_pi_series(acb_poly_t c, const acb_poly_t h, long n, long prec) + + Compute the respective trigonometric functions of the input + multiplied by `\pi`. + .. function:: void _acb_poly_tan_series(acb_ptr g, acb_srcptr h, long hlen, long len, long prec) .. function:: void acb_poly_tan_series(acb_poly_t g, const acb_poly_t h, long n, long prec) diff --git a/doc/source/arb_poly.rst b/doc/source/arb_poly.rst index 27dc2b5f..515d1782 100644 --- a/doc/source/arb_poly.rst +++ b/doc/source/arb_poly.rst @@ -766,13 +766,13 @@ Powers and elementary functions The underscore methods support aliasing and allow the input to be shorter than the output, but require the lengths to be nonzero. -.. function:: void _arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, long n, long prec) +.. function:: void _arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, long n, long prec, int times_pi) -.. function:: void arb_poly_sin_cos_series_basecase(arb_poly_t s, arb_poly_t c, const arb_poly_t h, long n, long prec) +.. function:: void arb_poly_sin_cos_series_basecase(arb_poly_t s, arb_poly_t c, const arb_poly_t h, long n, long prec, int times_pi) -.. function:: void _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, long n, long prec) +.. function:: void _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, long n, long prec, int times_pi) -.. function:: void arb_poly_sin_cos_series_tangent(arb_poly_t s, arb_poly_t c, const arb_poly_t h, long n, long prec) +.. function:: void arb_poly_sin_cos_series_tangent(arb_poly_t s, arb_poly_t c, const arb_poly_t h, long n, long prec, int times_pi) .. function:: void _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, long n, long prec) @@ -797,6 +797,9 @@ Powers and elementary functions The default version automatically selects between the *basecase* and *tangent* algorithms depending on the input. + The *basecase* and *tangent* versions take a flag *times_pi* + specifying that the input is to be multiplied by `\pi`. + The underscore methods support aliasing and require the lengths to be nonzero. .. function:: void _arb_poly_sin_series(arb_ptr s, arb_srcptr h, long hlen, long n, long prec) @@ -811,6 +814,21 @@ Powers and elementary functions simply wrap :func:`_arb_poly_sin_cos_series`. The underscore methods support aliasing and require the lengths to be nonzero. +.. function:: void _arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, arb_srcptr h, long hlen, long n, long prec) + +.. function:: void arb_poly_sin_cos_pi_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, long n, long prec) + +.. function:: void _arb_poly_sin_pi_series(arb_ptr s, arb_srcptr h, long hlen, long n, long prec) + +.. function:: void arb_poly_sin_pi_series(arb_poly_t s, const arb_poly_t h, long n, long prec) + +.. function:: void _arb_poly_cos_pi_series(arb_ptr c, arb_srcptr h, long hlen, long n, long prec) + +.. function:: void arb_poly_cos_pi_series(arb_poly_t c, const arb_poly_t h, long n, long prec) + + Compute the respective trigonometric functions of the input + multiplied by `\pi`. + .. function:: void _arb_poly_tan_series(arb_ptr g, arb_srcptr h, long hlen, long len, long prec) .. function:: void arb_poly_tan_series(arb_poly_t g, const arb_poly_t h, long n, long prec)