Skip to content

Commit d7f45cb

Browse files
fredrik-johanssonfredrik-johansson
authored
Merge pull request #2139 from fredrik-johansson/sinc
sinc_series: allow intervals containing 0
2 parents 1f9df9e + a6ff68c commit d7f45cb

File tree

3 files changed

+196
-1
lines changed

3 files changed

+196
-1
lines changed

‎src/acb_poly/sinc_series.c

+118-1
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
/*
2-
Copyright (C) 2016 Fredrik Johansson
2+
Copyright (C) 2016, 2025 Fredrik Johansson
33
44
This file is part of FLINT.
55
@@ -9,8 +9,118 @@
99
(at your option) any later version. See <https://www.gnu.org/licenses/>.
1010
*/
1111

12+
#include "mag.h"
1213
#include "acb_poly.h"
1314

15+
void
16+
_acb_sinc_jet_zero(acb_ptr res, const acb_t z, slong len, slong prec)
17+
{
18+
mag_ptr D;
19+
slong n;
20+
int is_real, is_imag;
21+
mag_t zm, err, fac;
22+
arb_t wide;
23+
24+
is_real = acb_is_real(z);
25+
is_imag = arb_is_zero(acb_realref(z));
26+
27+
mag_init(zm);
28+
mag_init(err);
29+
mag_init(fac);
30+
arb_init(wide);
31+
32+
acb_get_mag(zm, z);
33+
34+
D = _mag_vec_init(len + 1);
35+
36+
/* Compute sequence of derivative bounds
37+
38+
sinc^{(n)}(z) = int_0^1 t^n cos(t z + n pi / 2) dt
39+
40+
|sinc^{(n)}(z)| <= cosh(im(z)) min(1/|z|, 1/(n+1))
41+
|sinc^{(n)}(z)| <= cosh(im(z)) min(1/|z|, 1/(n+1), |z|/(n+2)) (odd n)
42+
*/
43+
{
44+
mag_t C, b1, b2, b3;
45+
46+
mag_init(C);
47+
mag_init(b1);
48+
mag_init(b2);
49+
mag_init(b3);
50+
51+
/* C = cosh(im(z)) */
52+
arb_get_mag(C, acb_imagref(z));
53+
mag_cosh(C, C);
54+
55+
/* b1 = 1/|z| */
56+
acb_get_mag_lower(b1, z);
57+
mag_inv(b1, b1);
58+
59+
for (n = 0; n <= len; n++)
60+
{
61+
/* b2 = 1/(n+1) */
62+
mag_one(b2);
63+
mag_div_ui(b2, b2, n + 1);
64+
mag_min(b2, b1, b2);
65+
66+
/* b3 = |z|/(n+2) */
67+
if (n % 2 == 1)
68+
{
69+
mag_div_ui(b3, zm, n + 2);
70+
mag_min(b2, b2, b3);
71+
}
72+
73+
mag_mul(D + n, C, b2);
74+
}
75+
76+
mag_clear(C);
77+
mag_clear(b1);
78+
mag_clear(b2);
79+
mag_clear(b3);
80+
}
81+
82+
mag_one(fac);
83+
84+
for (n = 0; n < len; n++)
85+
{
86+
/* sinc^{(n)}(0) / n! */
87+
if (n == 0)
88+
acb_one(res + n);
89+
else if (n % 2 == 1)
90+
acb_zero(res + n);
91+
else
92+
acb_div_si(res + n, res + n - 2, -n * (n + 1), prec);
93+
94+
if (n > 1)
95+
mag_div_ui(fac, fac, n);
96+
97+
/* |sinc^{(n)}(0 + eps) - sinc^{(n)}(0)| / n! <= |eps| * |sinc^{(n+1)}(z)| / n!, z = (+/- eps) */
98+
mag_mul(err, zm, D + n + 1);
99+
mag_mul(err, err, fac);
100+
acb_add_error_mag(res + n, err);
101+
102+
/* sinc^{(n)}(+/- eps) is a possibly better enclosure */
103+
arb_zero(wide);
104+
mag_mul(err, D + n, fac);
105+
arb_add_error_mag(wide, err);
106+
arb_intersection(acb_realref(res + n), acb_realref(res + n), wide, prec);
107+
arb_intersection(acb_imagref(res + n), acb_imagref(res + n), wide, prec);
108+
109+
if (is_real || (is_imag && (n % 2 == 0)))
110+
arb_zero(acb_imagref(res + n));
111+
112+
if (is_imag && (n % 2 == 1))
113+
arb_zero(acb_realref(res + n));
114+
}
115+
116+
_mag_vec_clear(D, len + 1);
117+
118+
mag_clear(zm);
119+
mag_clear(err);
120+
mag_clear(fac);
121+
arb_clear(wide);
122+
}
123+
14124
void
15125
_acb_poly_sinc_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec)
16126
{
@@ -35,6 +145,13 @@ _acb_poly_sinc_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec)
35145
_acb_poly_sin_series(t, u, hlen, n + 1, prec);
36146
_acb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec);
37147
}
148+
else if (acb_contains_zero(h))
149+
{
150+
_acb_sinc_jet_zero(t, h, n, prec);
151+
/* compose with nonconstant part */
152+
acb_zero(u);
153+
_acb_poly_compose_series(g, t, n, u, hlen, n, prec);
154+
}
38155
else
39156
{
40157
_acb_poly_sin_series(t, u, hlen, n, prec);

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

+54
Original file line numberDiff line numberDiff line change
@@ -83,5 +83,59 @@ TEST_FUNCTION_START(acb_poly_sinc_series, state)
8383
acb_poly_clear(d);
8484
}
8585

86+
/* check intervals containing 0 */
87+
for (iter = 0; iter < 100 * flint_test_multiplier(); iter++)
88+
{
89+
acb_t z1, z2;
90+
acb_poly_t f1, f2, g1, g2;
91+
slong len, prec;
92+
93+
acb_init(z1);
94+
acb_init(z2);
95+
96+
acb_poly_init(f1);
97+
acb_poly_init(f2);
98+
acb_poly_init(g1);
99+
acb_poly_init(g2);
100+
101+
prec = 2 + n_randint(state, 300);
102+
len = 1 + n_randint(state, 5);
103+
104+
acb_randtest(z1, state, 300, 5);
105+
acb_mul_2exp_si(z1, z1, -(slong) n_randint(state, 100));
106+
if (n_randint(state, 2))
107+
arb_zero(acb_realref(z1));
108+
109+
arb_get_mag(arb_radref(acb_realref(z2)), acb_realref(z1));
110+
arb_get_mag(arb_radref(acb_imagref(z2)), acb_imagref(z1));
111+
112+
acb_poly_set_coeff_si(f1, 1, n_randint(state, 2) ? 1 : -1);
113+
acb_poly_set(f2, f1);
114+
acb_poly_set_coeff_acb(f1, 0, z1);
115+
acb_poly_set_coeff_acb(f2, 0, z2);
116+
117+
acb_poly_sinc_series(g1, f1, len, prec);
118+
acb_poly_sinc_series(g2, f2, len, prec);
119+
120+
if (!acb_poly_overlaps(g1, g2))
121+
{
122+
flint_printf("FAIL: overlap (0)\n\n");
123+
flint_printf("iter = %wd\n", iter);
124+
flint_printf("f1 = "); acb_poly_printd(f1, prec / 3.33); flint_printf("\n\n");
125+
flint_printf("f2 = "); acb_poly_printd(f2, prec / 3.33); flint_printf("\n\n");
126+
flint_printf("g1 = "); acb_poly_printd(g1, prec / 3.33); flint_printf("\n\n");
127+
flint_printf("g2 = "); acb_poly_printd(g2, prec / 3.33); flint_printf("\n\n");
128+
flint_abort();
129+
}
130+
131+
acb_clear(z1);
132+
acb_clear(z2);
133+
134+
acb_poly_clear(f1);
135+
acb_poly_clear(f2);
136+
acb_poly_clear(g1);
137+
acb_poly_clear(g2);
138+
}
139+
86140
TEST_FUNCTION_END(state);
87141
}

‎src/arb_poly/sinc_series.c

+24
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,23 @@
1010
*/
1111

1212
#include "arb_poly.h"
13+
#include "acb_poly.h"
14+
15+
void _acb_sinc_jet_zero(acb_ptr res, const acb_t z, slong len, slong prec);
16+
17+
void
18+
_arb_sinc_jet_zero(arb_ptr res, const arb_t z, slong len, slong prec)
19+
{
20+
acb_ptr t;
21+
slong i;
22+
t = _acb_vec_init(len + 1);
23+
acb_set_arb(t + len, z);
24+
_acb_sinc_jet_zero(t, t + len, len, prec);
25+
for (i = 0; i < len; i++)
26+
arb_swap(res + i, acb_realref(t + i));
27+
_acb_vec_clear(t, len + 1);
28+
29+
}
1330

1431
void
1532
_arb_poly_sinc_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec)
@@ -35,6 +52,13 @@ _arb_poly_sinc_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec)
3552
_arb_poly_sin_series(t, u, hlen, n + 1, prec);
3653
_arb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec);
3754
}
55+
else if (arb_contains_zero(h))
56+
{
57+
_arb_sinc_jet_zero(t, h, n, prec);
58+
/* compose with nonconstant part */
59+
arb_zero(u);
60+
_arb_poly_compose_series(g, t, n, u, hlen, n, prec);
61+
}
3862
else
3963
{
4064
_arb_poly_sin_series(t, u, hlen, n, prec);

0 commit comments

Comments
 (0)