1 /*
2 Copyright (C) 2017 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 /* Computes sin(x) or cos(x) using Taylor series truncated at x^N exclusive.
15 Computes error bound automatically. Does not allow aliasing of s and x. */
16 void
arb_sin_cos_taylor_sum_rs(arb_t s,const arb_t x,slong N,int cosine,slong prec)17 arb_sin_cos_taylor_sum_rs(arb_t s, const arb_t x, slong N, int cosine, slong prec)
18 {
19 mag_t err;
20 mag_init(err);
21
22 arb_get_mag(err, x);
23 mag_exp_tail(err, err, N);
24
25 if (N == 0 || (!cosine && N == 1))
26 {
27 arb_zero(s);
28 }
29 else if (cosine && N <= 2)
30 {
31 arb_one(s);
32 }
33 else if (!cosine && N <= 3) /* x */
34 {
35 arb_set_round(s, x, prec);
36 }
37 else if (cosine && N <= 4) /* 1 - x^2/2 */
38 {
39 arb_mul(s, x, x, prec / 2 + 4);
40 arb_mul_2exp_si(s, s, -1);
41 arb_sub_ui(s, s, 1, prec);
42 arb_neg(s, s);
43 }
44 else if (!cosine && N <= 5) /* x - x^3/6 */
45 {
46 arb_mul(s, x, x, prec / 2 + 4);
47 arb_div_ui(s, s, 6, prec / 2 + 4);
48 arb_mul(s, s, x, prec / 2 + 4);
49 arb_sub(s, x, s, prec);
50 }
51 else
52 {
53 arb_ptr tpow;
54 slong j, k, m, M, tp, xmag;
55 mp_limb_t c, d, chi, clo;
56
57 xmag = arf_abs_bound_lt_2exp_si(arb_midref(x));
58
59 /* Convert to order as a series in x^2. */
60 if (cosine)
61 M = (N + 1) / 2;
62 else
63 M = N / 2;
64
65 m = n_sqrt(M);
66
67 /* not intended (and not 32-bit safe...) */
68 if (M > 30000)
69 {
70 flint_abort();
71 }
72
73 tpow = _arb_vec_init(m + 1);
74
75 arb_mul(s, x, x, prec);
76 _arb_vec_set_powers(tpow, s, m + 1, prec);
77 arb_zero(s);
78
79 c = 1;
80 j = (M - 1) % m;
81
82 for (k = M - 1; k >= 0; k--)
83 {
84 tp = prec - 2 * k * (-xmag) + 10;
85 tp = FLINT_MAX(tp, 2);
86 tp = FLINT_MIN(tp, prec);
87
88 if (cosine)
89 d = (2 * k - 1) * (2 * k);
90 else
91 d = (2 * k) * (2 * k + 1);
92
93 if (k != 0)
94 {
95 umul_ppmm(chi, clo, c, d);
96
97 if (chi != 0)
98 {
99 arb_div_ui(s, s, c, tp);
100 c = 1;
101 }
102 }
103
104 if (k % 2 == 0)
105 arb_addmul_ui(s, tpow + j, c, tp);
106 else
107 arb_submul_ui(s, tpow + j, c, tp);
108
109 if (k != 0)
110 {
111 c *= d;
112
113 if (j == 0)
114 {
115 arb_mul(s, s, tpow + m, tp);
116 j = m - 1;
117 }
118 else
119 {
120 j--;
121 }
122 }
123 }
124
125 arb_div_ui(s, s, c, prec);
126 if (!cosine)
127 arb_mul(s, s, x, prec);
128
129 _arb_vec_clear(tpow, m + 1);
130 }
131
132 arb_add_error_mag(s, err);
133 mag_clear(err);
134 }
135
136 void
arb_sin_cos_arf_rs_generic(arb_t res_sin,arb_t res_cos,const arf_t x,slong prec)137 arb_sin_cos_arf_rs_generic(arb_t res_sin, arb_t res_cos, const arf_t x, slong prec)
138 {
139 slong q, xmag, wp, k, N;
140 arb_t s, t;
141 int negate;
142
143 if (arf_is_zero(x))
144 {
145 if (res_sin != NULL)
146 arb_zero(res_sin);
147 if (res_cos != NULL)
148 arb_one(res_cos);
149 return;
150 }
151
152 xmag = arf_abs_bound_lt_2exp_si(x);
153
154 /* x + O(x^3), 1 + O(x^2) */
155 if (xmag < -(prec / 2) - 4)
156 {
157 arb_init(t);
158 arf_set(arb_midref(t), x);
159 if (res_sin != NULL)
160 arb_sin_cos_taylor_sum_rs(res_sin, t, 3, 0, prec);
161 if (res_cos != NULL)
162 arb_sin_cos_taylor_sum_rs(res_cos, t, 2, 1, prec);
163 arb_clear(t);
164 return;
165 }
166
167 xmag = FLINT_MAX(xmag, -prec);
168
169 /* could include sanity test, but we assume that the function
170 gets called with a valid |x| < pi/2
171 if (xmag > 0)
172 {
173 if (res_sin != NULL)
174 arb_indeterminate(res_sin);
175 if (res_cos != NULL)
176 arb_indeterminate(res_cos);
177 return;
178 }
179 */
180
181 arb_init(s);
182 arb_init(t);
183
184 negate = arf_sgn(x) < 0;
185
186 /* generic tuning value */
187 q = 4.5 * pow(prec, 0.2);
188 q = FLINT_MAX(q, 6);
189 /* adjust to magnitude */
190 q = FLINT_MAX(0, xmag + q);
191 /* don't do a redundant square root */
192 if (q <= 2)
193 q = 0;
194
195 wp = prec + 10 + 2 * FLINT_BIT_COUNT(prec);
196
197 /* t = x/2^q */
198 arf_mul_2exp_si(arb_midref(t), x, -q);
199
200 if (q == 0 && res_sin != NULL)
201 {
202 /* compute cos from sin since the square root has less cancellation */
203 wp += (-xmag);
204 N = _arb_exp_taylor_bound(xmag, wp);
205 arb_sin_cos_taylor_sum_rs(s, t, N, 0, wp);
206
207 if (res_sin != NULL)
208 arb_set_round(res_sin, s, prec);
209
210 if (res_cos != NULL)
211 {
212 arb_mul(t, s, s, wp);
213 arb_sub_ui(t, t, 1, wp);
214 arb_neg(t, t);
215 arb_sqrt(res_cos, t, prec);
216 }
217 }
218 else
219 {
220 /* compute sin from cos */
221 wp = prec + 10 + 2 * FLINT_BIT_COUNT(prec);
222 wp += 2 * (q - xmag); /* todo: too much when only computing cos? */
223
224 N = _arb_exp_taylor_bound(xmag - q, wp);
225
226 arb_sin_cos_taylor_sum_rs(s, t, N, 1, wp);
227
228 for (k = 0; k < q; k++)
229 {
230 arb_mul(s, s, s, wp);
231 arb_mul_2exp_si(s, s, 1);
232 arb_sub_ui(s, s, 1, wp);
233 }
234
235 if (res_cos != NULL)
236 arb_set_round(res_cos, s, prec);
237
238 if (res_sin != NULL)
239 {
240 arb_mul(s, s, s, wp);
241 arb_sub_ui(s, s, 1, wp);
242 arb_neg(s, s);
243 arb_sqrtpos(res_sin, s, prec);
244 if (negate)
245 arb_neg(res_sin, res_sin);
246 }
247 }
248
249 arb_clear(s);
250 arb_clear(t);
251 }
252
253 void
arb_sin_cos_arf_generic(arb_t res_sin,arb_t res_cos,const arf_t x,slong prec)254 arb_sin_cos_arf_generic(arb_t res_sin, arb_t res_cos, const arf_t x, slong prec)
255 {
256 arb_t pi4, t, u, v;
257 fmpz_t q;
258 slong wp, mag;
259 int octant, swapsincos, sinnegative, cosnegative, negative;
260
261 mag = arf_abs_bound_lt_2exp_si(x);
262
263 if (mag > FLINT_MAX(65536, 4 * prec))
264 {
265 if (res_sin != NULL)
266 arb_zero_pm_one(res_sin);
267 if (res_cos != NULL)
268 arb_zero_pm_one(res_cos);
269 }
270 else if (mag <= 0) /* todo: compare with pi/4-eps instead? */
271 {
272 if (prec < 90000 || mag < -prec / 16 ||
273 /* rs is faster for even smaller prec/N than this but has high memory usage */
274 (prec < 100000000 && mag < -prec / 128))
275 arb_sin_cos_arf_rs_generic(res_sin, res_cos, x, prec);
276 else
277 arb_sin_cos_arf_bb(res_sin, res_cos, x, prec);
278 }
279 else
280 {
281 arb_init(pi4);
282 arb_init(t);
283 arb_init(u);
284 arb_init(v);
285 fmpz_init(q);
286
287 wp = prec + mag + 10;
288 negative = arf_sgn(x) < 0;
289
290 arb_const_pi(pi4, wp);
291 arb_mul_2exp_si(pi4, pi4, -2);
292 arb_set_arf(t, x);
293 arb_abs(t, t);
294
295 arb_set_round(v, t, mag + 10);
296 arb_set_round(u, pi4, mag + 10);
297 arb_div(u, v, u, mag + 10);
298
299 arf_get_fmpz(q, arb_midref(u), ARF_RND_DOWN);
300 arb_submul_fmpz(t, pi4, q, wp);
301 octant = fmpz_fdiv_ui(q, 8);
302 if (octant & 1)
303 arb_sub(t, pi4, t, wp);
304
305 arb_clear(pi4);
306 arb_clear(u);
307 arb_clear(v);
308
309 sinnegative = (octant >= 4) ^ negative;
310 cosnegative = (octant >= 2 && octant <= 5);
311 swapsincos = (octant == 1 || octant == 2 || octant == 5 || octant == 6);
312
313 /* guard against infinite recursion */
314 if (arf_cmpabs_2exp_si(arb_midref(t), 0) > 0)
315 {
316 flint_printf("mod pi/4 reduction unexpectedly failed!\n");
317 flint_abort();
318 }
319
320 /* todo: allow NULL in arb_sin_cos and simplify here */
321 if (swapsincos)
322 {
323 if (res_sin != NULL && res_cos != NULL)
324 arb_sin_cos(res_cos, res_sin, t, prec);
325 else if (res_sin != NULL)
326 arb_cos(res_sin, t, prec);
327 else
328 arb_sin(res_cos, t, prec);
329 }
330 else
331 {
332 if (res_sin != NULL && res_cos != NULL)
333 arb_sin_cos(res_sin, res_cos, t, prec);
334 else if (res_sin != NULL)
335 arb_sin(res_sin, t, prec);
336 else
337 arb_cos(res_cos, t, prec);
338 }
339
340 if (sinnegative && res_sin != NULL)
341 arb_neg(res_sin, res_sin);
342 if (cosnegative && res_cos != NULL)
343 arb_neg(res_cos, res_cos);
344
345 arb_clear(t);
346 fmpz_clear(q);
347 }
348 }
349
350