1 /*
2 Copyright (C) 2018 Fredrik Johansson
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 #define ONE_OVER_PI 0.31830988618379067154
15 #define PI 3.1415926535897932385
16
17 /* We use doubles in a way that keeps the final error <= EPSILON. */
18 #define EPSILON ldexp(1.0, -30)
19
20 /* For |x| <= 2^MAX_EXP, doubles can be used directly for quadrant
21 reduction. Note: 2^MAX_EXP must also fit in an int. */
22 #define MAX_EXP 20
23
24 /* Lookup tables, steps of 1/16. */
25 static const double sin_tab[] = {
26 0.0,0.062459317842380198585,0.12467473338522768996,
27 0.18640329676226988455,0.24740395925452292960,0.30743851458038085067,
28 0.36627252908604756137,0.42367625720393801036,0.47942553860420300027,
29 0.53330267353602017333,0.58509727294046215481,0.63460708001526929685,
30 0.68163876002333416673,0.72600865526071254966,0.76754350223602703963,
31 0.80608110826069299518,0.84147098480789650665,0.87357493516707112023,
32 0.90226759409909516292,0.92743691738486767172,0.94898461935558621435,
33 0.96682655669618022997,0.98089305702315569609,0.99112919095376166988,
34 0.99749498660405443094,0.99996558567824887530,
35 };
36
37 static const double cos_tab[] = {
38 1.0,0.99804751070009914963,0.99219766722932905315,
39 0.98247331310125525749,0.96891242171064478414,0.95156794804817220215,
40 0.93050762191231429115,0.90581368342593642074,0.87758256189037271612,
41 0.84592449923106795446,0.81096311950521790219,0.77283494615247154481,
42 0.73168886887382088631,0.68768556222050484451,0.64099685816332513036,
43 0.59180507509247750546,0.54030230586813971740,0.48668966770196333087,
44 0.43117651679866617655,0.37397963082453319229,0.31532236239526866545,
45 0.25543376688881169791,0.19454770798898718445,0.13290194445282520566,
46 0.070737201667702910088,0.0082962316238583774779,
47 };
48
49 /* Computes sin(a) and cos(a) and also sets q to the (approximate)
50 quadrant of a. The absolute error of the straightforward algorithm
51 below can be bounded as the error in the mod pi/2 reduction plus
52 the negligible contribution (O(1) * 2^-53) of all other steps.
53 Since |a| < 2^20 by assumption, we have 3 bits of margin for the
54 final error bound of 2^-30 for the entire algorithm. */
55 static void
sin_cos(double * sin_a,double * cos_a,int * q,double a)56 sin_cos(double * sin_a, double * cos_a, int * q, double a)
57 {
58 double as, ac, t, b, bs, bc, v;
59 int i, qa;
60
61 *q = qa = floor(a * (2.0 * ONE_OVER_PI));
62 a = a - qa * (0.5 * PI);
63
64 if (a < 0.0)
65 a = 0.0;
66 if (a > 0.5 * PI)
67 a = 0.5 * PI;
68
69 i = a * 16.0;
70
71 if (i < 0 || i > 25)
72 flint_abort();
73
74 as = sin_tab[i];
75 ac = cos_tab[i];
76
77 b = a - i * (1 / 16.0);
78 v = b * b;
79
80 /* Just use the Taylor series. */
81 bs = 2.7557319223985890653e-6 * v;
82 bs = (-0.0001984126984126984127 + bs) * v;
83 bs = (0.0083333333333333333333 + bs) * v;
84 bs = (-0.16666666666666666667 + bs) * v;
85 bs = (1.0 + bs) * b;
86
87 bc = -2.7557319223985890653e-7 * v;
88 bc = (0.000024801587301587301587 + bc) * v;
89 bc = (-0.0013888888888888888889 + bc) * v;
90 bc = (0.041666666666666666667 + bc) * v;
91 bc = (-0.5 + bc) * v;
92 bc = 1.0 + bc;
93
94 t = as * bc + ac * bs;
95 ac = ac * bc - as * bs;
96 as = t;
97
98 if ((qa & 3) == 0)
99 {
100 *sin_a = as;
101 *cos_a = ac;
102 }
103 else if ((qa & 3) == 1)
104 {
105 *sin_a = ac;
106 *cos_a = -as;
107 }
108 else if ((qa & 3) == 2)
109 {
110 *sin_a = -as;
111 *cos_a = -ac;
112 }
113 else
114 {
115 *sin_a = -ac;
116 *cos_a = as;
117 }
118 }
119
120 /* FIXME: this is the bottleneck -- an mpn version would be better */
121 static void
_arb_mod_2pi(arb_t x,slong mag)122 _arb_mod_2pi(arb_t x, slong mag)
123 {
124 arb_t t;
125 arf_t q;
126
127 arf_init(q);
128 arb_init(t);
129
130 arb_const_pi(t, mag + 53);
131 arb_mul_2exp_si(t, t, 1);
132 arf_div(q, arb_midref(x), arb_midref(t), mag + 10, ARF_RND_NEAR);
133 arf_floor(q, q);
134 arb_submul_arf(x, t, q, 53);
135
136 arf_clear(q);
137 arb_clear(t);
138 }
139
140 void
_arb_sin_cos_wide(arb_t sinx,arb_t cosx,const arf_t xmid,const mag_t xrad,slong prec)141 _arb_sin_cos_wide(arb_t sinx, arb_t cosx, const arf_t xmid, const mag_t xrad, slong prec)
142 {
143 double m, a, b, r, cos_min, cos_max, sin_min, sin_max;
144 double as, ac, bs, bc;
145 int i, qa, qb;
146 slong mag;
147
148 mag = arf_abs_bound_lt_2exp_si(xmid);
149
150 if (mag > FLINT_MAX(65536, 4 * prec) || mag_cmp_2exp_si(xrad, 3) >= 0)
151 {
152 if (sinx != NULL)
153 arb_zero_pm_one(sinx);
154 if (cosx != NULL)
155 arb_zero_pm_one(cosx);
156 return;
157 }
158 else if (mag <= MAX_EXP)
159 {
160 m = arf_get_d(xmid, ARF_RND_DOWN);
161 r = mag_get_d(xrad);
162 }
163 else
164 {
165 arb_t t;
166 arb_init(t);
167 arf_set(arb_midref(t), xmid);
168 mag_set(arb_radref(t), xrad);
169 _arb_mod_2pi(t, mag);
170
171 /* this should not happen */
172 if (arf_cmpabs_2exp_si(arb_midref(t), 5) > 0 ||
173 mag_cmp_2exp_si(arb_radref(t), 5) > 0)
174 {
175 flint_printf("unexpected precision loss in sin_cos_wide\n");
176
177 if (sinx != NULL)
178 arb_zero_pm_one(sinx);
179 if (cosx != NULL)
180 arb_zero_pm_one(cosx);
181 arb_clear(t);
182 return;
183 }
184
185 m = arf_get_d(arb_midref(t), ARF_RND_DOWN);
186 r = mag_get_d(arb_radref(t));
187
188 arb_clear(t);
189 }
190
191 a = m - r;
192 b = m + r;
193
194 sin_cos(&as, &ac, &qa, a);
195 sin_cos(&bs, &bc, &qb, b);
196
197 sin_min = FLINT_MIN(as, bs);
198 sin_max = FLINT_MAX(as, bs);
199 cos_min = FLINT_MIN(ac, bc);
200 cos_max = FLINT_MAX(ac, bc);
201
202 /* Handle the quadrant crossings. */
203 for (i = qa; i < qb; i++)
204 {
205 if ((i & 3) == 1) cos_min = -1.0;
206 if ((i & 3) == 3) cos_max = 1.0;
207 if ((i & 3) == 2) sin_min = -1.0;
208 if ((i & 3) == 0) sin_max = 1.0;
209 }
210
211 if (sinx != NULL)
212 {
213 a = (sin_max + sin_min) * 0.5;
214 r = (sin_max - sin_min) * 0.5 + EPSILON;
215
216 arf_set_d(arb_midref(sinx), a);
217 mag_set_d(arb_radref(sinx), r);
218 arb_set_round(sinx, sinx, prec);
219 }
220
221 if (cosx != NULL)
222 {
223 a = (cos_max + cos_min) * 0.5;
224 r = (cos_max - cos_min) * 0.5 + EPSILON;
225
226 arf_set_d(arb_midref(cosx), a);
227 mag_set_d(arb_radref(cosx), r);
228 arb_set_round(cosx, cosx, prec);
229 }
230 }
231
232 void
arb_sin_cos_wide(arb_t sinx,arb_t cosx,const arb_t x,slong prec)233 arb_sin_cos_wide(arb_t sinx, arb_t cosx, const arb_t x, slong prec)
234 {
235 _arb_sin_cos_wide(sinx, cosx, arb_midref(x), arb_radref(x), prec);
236 }
237