Skip to content

Commit f0b28ad

Browse files
fredrik-johanssonfredrik-johansson
authored
Merge pull request #2142 from fredrik-johansson/sinc2
Handle 0-containing constant term in arb_poly_sinc_pi_series; add acb_poly_sinc_pi_series
2 parents d7f45cb + 3a4131b commit f0b28ad

File tree

6 files changed

+198
-1
lines changed

6 files changed

+198
-1
lines changed

‎doc/source/acb_poly.rst

+7
Original file line numberDiff line numberDiff line change
@@ -848,6 +848,13 @@ Elementary functions
848848
Sets *s* to the sinc function of the power series *h*, truncated
849849
to length *n*.
850850

851+
.. function:: void _acb_poly_sinc_pi_series(acb_ptr s, acb_srcptr h, slong hlen, slong n, slong prec)
852+
853+
.. function:: void acb_poly_sinc_pi_series(acb_poly_t s, const acb_poly_t h, slong n, slong prec)
854+
855+
Compute the sinc function of the input multiplied by `\pi`.
856+
857+
851858
Lambert W function
852859
-------------------------------------------------------------------------------
853860

‎src/acb_poly.h

+2
Original file line numberDiff line numberDiff line change
@@ -564,6 +564,8 @@ void acb_poly_tan_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec);
564564

565565
void _acb_poly_sinc_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec);
566566
void acb_poly_sinc_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec);
567+
void _acb_poly_sinc_pi_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec);
568+
void acb_poly_sinc_pi_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec);
567569

568570
void _acb_poly_lambertw_series(acb_ptr res, acb_srcptr z, slong zlen, const fmpz_t k, int flags, slong len, slong prec);
569571
void acb_poly_lambertw_series(acb_poly_t res, const acb_poly_t z, const fmpz_t k, int flags, slong len, slong prec);

‎src/acb_poly/sinc_pi_series.c

+85
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
/*
2+
Copyright (C) 2016 Fredrik Johansson
3+
Copyright (C) 2019 D.H.J. Polymath
4+
5+
This file is part of FLINT.
6+
7+
FLINT is free software: you can redistribute it and/or modify it under
8+
the terms of the GNU Lesser General Public License (LGPL) as published
9+
by the Free Software Foundation; either version 3 of the License, or
10+
(at your option) any later version. See <https://www.gnu.org/licenses/>.
11+
*/
12+
13+
#include "acb_poly.h"
14+
15+
void
16+
_acb_poly_sinc_pi_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec)
17+
{
18+
hlen = FLINT_MIN(hlen, n);
19+
20+
if (hlen == 1)
21+
{
22+
acb_sinc_pi(g, h, prec);
23+
_acb_vec_zero(g + 1, n - 1);
24+
}
25+
else
26+
{
27+
acb_t pi;
28+
acb_ptr t, u;
29+
30+
acb_init(pi);
31+
t = _acb_vec_init(n + 1);
32+
u = _acb_vec_init(hlen);
33+
34+
acb_const_pi(pi, prec);
35+
_acb_vec_set(u, h, hlen);
36+
37+
if (acb_is_zero(h))
38+
{
39+
_acb_poly_sin_pi_series(t, u, hlen, n + 1, prec);
40+
_acb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec);
41+
_acb_vec_scalar_div(g, g, n, pi, prec);
42+
}
43+
else if (acb_contains_zero(h))
44+
{
45+
_acb_vec_scalar_mul(u, h, hlen, pi, prec);
46+
_acb_poly_sinc_series(g, u, hlen, n, prec);
47+
}
48+
else
49+
{
50+
_acb_poly_sin_pi_series(t, u, hlen, n, prec);
51+
_acb_poly_div_series(g, t, n, u, hlen, n, prec);
52+
_acb_vec_scalar_div(g, g, n, pi, prec);
53+
}
54+
55+
acb_clear(pi);
56+
_acb_vec_clear(t, n + 1);
57+
_acb_vec_clear(u, hlen);
58+
}
59+
}
60+
61+
void
62+
acb_poly_sinc_pi_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec)
63+
{
64+
slong hlen = h->length;
65+
66+
if (n == 0)
67+
{
68+
acb_poly_zero(g);
69+
return;
70+
}
71+
72+
if (hlen == 0)
73+
{
74+
acb_poly_one(g);
75+
return;
76+
}
77+
78+
if (hlen == 1)
79+
n = 1;
80+
81+
acb_poly_fit_length(g, n);
82+
_acb_poly_sinc_pi_series(g->coeffs, h->coeffs, hlen, n, prec);
83+
_acb_poly_set_length(g, n);
84+
_acb_poly_normalise(g);
85+
}

‎src/acb_poly/test/main.c

+2
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@
7777
#include "t-sin_cos_pi_series.c"
7878
#include "t-sin_cos_series.c"
7979
#include "t-sinc_series.c"
80+
#include "t-sinc_pi_series.c"
8081
#include "t-sinh_cosh_series.c"
8182
#include "t-sin_pi_series.c"
8283
#include "t-sin_series_cos_series.c"
@@ -159,6 +160,7 @@ test_struct tests[] =
159160
TEST_FUNCTION(acb_poly_sin_cos_pi_series),
160161
TEST_FUNCTION(acb_poly_sin_cos_series),
161162
TEST_FUNCTION(acb_poly_sinc_series),
163+
TEST_FUNCTION(acb_poly_sinc_pi_series),
162164
TEST_FUNCTION(acb_poly_sinh_cosh_series),
163165
TEST_FUNCTION(acb_poly_sin_pi_series),
164166
TEST_FUNCTION(acb_poly_sin_series_cos_series),

‎src/acb_poly/test/t-sinc_pi_series.c

+95
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
/*
2+
Copyright (C) 2012, 2013 Fredrik Johansson
3+
Copyright (C) 2019 D.H.J. Polymath
4+
5+
This file is part of FLINT.
6+
7+
FLINT is free software: you can redistribute it and/or modify it under
8+
the terms of the GNU Lesser General Public License (LGPL) as published
9+
by the Free Software Foundation; either version 3 of the License, or
10+
(at your option) any later version. See <https://www.gnu.org/licenses/>.
11+
*/
12+
13+
#include "test_helpers.h"
14+
#include "acb_poly.h"
15+
16+
TEST_FUNCTION_START(acb_poly_sinc_pi_series, state)
17+
{
18+
slong iter;
19+
20+
for (iter = 0; iter < 200 * 0.1 * flint_test_multiplier(); iter++)
21+
{
22+
slong m, n1, n2, rbits1, rbits2, rbits3, rbits4;
23+
acb_poly_t a, b, c, d;
24+
acb_t pi;
25+
26+
rbits1 = 2 + n_randint(state, 300);
27+
rbits2 = 2 + n_randint(state, 300);
28+
rbits3 = 2 + n_randint(state, 300);
29+
rbits4 = 2 + n_randint(state, 300);
30+
31+
m = n_randint(state, 15);
32+
n1 = n_randint(state, 15);
33+
n2 = n_randint(state, 15);
34+
35+
acb_poly_init(a);
36+
acb_poly_init(b);
37+
acb_poly_init(c);
38+
acb_poly_init(d);
39+
acb_init(pi);
40+
41+
acb_poly_randtest(a, state, m, rbits1, 10);
42+
acb_poly_randtest(b, state, 10, rbits1, 10);
43+
acb_poly_randtest(c, state, 10, rbits1, 10);
44+
45+
acb_poly_sinc_pi_series(b, a, n1, rbits2);
46+
acb_poly_sinc_pi_series(c, a, n2, rbits3);
47+
48+
acb_poly_set(d, b);
49+
acb_poly_truncate(d, FLINT_MIN(n1, n2));
50+
acb_poly_truncate(c, FLINT_MIN(n1, n2));
51+
52+
acb_const_pi(pi, rbits4);
53+
54+
if (!acb_poly_overlaps(c, d))
55+
{
56+
flint_printf("FAIL\n\n");
57+
flint_printf("n1 = %wd, n2 = %wd, bits2 = %wd, bits3 = %wd, bits4 = %wd\n", n1, n2, rbits2, rbits3, rbits4);
58+
flint_printf("a = "); acb_poly_printd(a, 50); flint_printf("\n\n");
59+
flint_printf("b = "); acb_poly_printd(b, 50); flint_printf("\n\n");
60+
flint_printf("c = "); acb_poly_printd(c, 50); flint_printf("\n\n");
61+
flint_abort();
62+
}
63+
64+
/* check pi x sinc_pi(x) = sin_pi(x) */
65+
acb_poly_mullow(c, b, a, n1, rbits2);
66+
acb_poly_scalar_mul(c, c, pi, rbits2);
67+
acb_poly_sin_pi_series(d, a, n1, rbits2);
68+
69+
if (!acb_poly_overlaps(c, d))
70+
{
71+
flint_printf("FAIL (functional equation)\n\n");
72+
flint_printf("a = "); acb_poly_printd(a, 15); flint_printf("\n\n");
73+
flint_printf("b = "); acb_poly_printd(b, 15); flint_printf("\n\n");
74+
flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
75+
flint_printf("d = "); acb_poly_printd(d, 15); flint_printf("\n\n");
76+
flint_abort();
77+
}
78+
79+
acb_poly_sinc_pi_series(a, a, n1, rbits2);
80+
81+
if (!acb_poly_overlaps(a, b))
82+
{
83+
flint_printf("FAIL (aliasing)\n\n");
84+
flint_abort();
85+
}
86+
87+
acb_poly_clear(a);
88+
acb_poly_clear(b);
89+
acb_poly_clear(c);
90+
acb_poly_clear(d);
91+
acb_clear(pi);
92+
}
93+
94+
TEST_FUNCTION_END(state);
95+
}

‎src/arb_poly/sinc_pi_series.c

+7-1
Original file line numberDiff line numberDiff line change
@@ -38,13 +38,19 @@ _arb_poly_sinc_pi_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong pre
3838
{
3939
_arb_poly_sin_pi_series(t, u, hlen, n + 1, prec);
4040
_arb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec);
41+
_arb_vec_scalar_div(g, g, n, pi, prec);
42+
}
43+
else if (arb_contains_zero(h))
44+
{
45+
_arb_vec_scalar_mul(u, h, hlen, pi, prec);
46+
_arb_poly_sinc_series(g, u, hlen, n, prec);
4147
}
4248
else
4349
{
4450
_arb_poly_sin_pi_series(t, u, hlen, n, prec);
4551
_arb_poly_div_series(g, t, n, u, hlen, n, prec);
52+
_arb_vec_scalar_div(g, g, n, pi, prec);
4653
}
47-
_arb_vec_scalar_div(g, g, n, pi, prec);
4854

4955
arb_clear(pi);
5056
_arb_vec_clear(t, n + 1);

0 commit comments

Comments
 (0)