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