1 /*
2 Copyright (C) 2015 Arb authors
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "arb.h"
13
14 void
_arb_sinc_derivative_bound(mag_t d,const arb_t x)15 _arb_sinc_derivative_bound(mag_t d, const arb_t x)
16 {
17 /* |f'(x)| < min(arb_get_mag(x), 1) / 2 */
18 mag_t r, one;
19 mag_init(r);
20 mag_init(one);
21 arb_get_mag(r, x);
22 mag_one(one);
23 mag_min(d, r, one);
24 mag_mul_2exp_si(d, d, -1);
25 mag_clear(r);
26 mag_clear(one);
27 }
28
29 void
_arb_sinc_direct(arb_t z,const arb_t x,slong prec)30 _arb_sinc_direct(arb_t z, const arb_t x, slong prec)
31 {
32 /* z = sin(x) / x */
33 slong wp;
34 arb_t y;
35 wp = prec + 2;
36 arb_init(y);
37 arb_sin(y, x, wp);
38 arb_div(z, y, x, prec);
39 arb_clear(y);
40 }
41
42 void
arb_sinc(arb_t z,const arb_t x,slong prec)43 arb_sinc(arb_t z, const arb_t x, slong prec)
44 {
45 mag_t c, r;
46 mag_init(c);
47 mag_init(r);
48 mag_set_ui_2exp_si(c, 5, -1);
49 arb_get_mag_lower(r, x);
50 if (mag_cmp(c, r) < 0)
51 {
52 /* x is not near the origin */
53 _arb_sinc_direct(z, x, prec);
54 }
55 else if (mag_cmp_2exp_si(arb_radref(x), 1) < 0)
56 {
57 /* determine error magnitude using the derivative bound */
58 if (arb_is_exact(x))
59 {
60 mag_zero(c);
61 }
62 else
63 {
64 _arb_sinc_derivative_bound(r, x);
65 mag_mul(c, arb_radref(x), r);
66 }
67
68 /* evaluate sinc at the midpoint of x */
69 if (arf_is_zero(arb_midref(x)))
70 {
71 arb_one(z);
72 }
73 else
74 {
75 arb_get_mid_arb(z, x);
76 _arb_sinc_direct(z, z, prec);
77 }
78
79 /* add the error */
80 mag_add(arb_radref(z), arb_radref(z), c);
81 }
82 else
83 {
84 /* x has a large radius and includes points near the origin */
85 arf_zero(arb_midref(z));
86 mag_one(arb_radref(z));
87 }
88
89 mag_clear(c);
90 mag_clear(r);
91 }
92