1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2020 Free Software Foundation, Inc.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <https://www.gnu.org/licenses/>.
18 */
19/****************************************************************/
20/* MODULE_NAME: sincos32.c */
21/* */
22/* FUNCTIONS: ss32 */
23/* cc32 */
24/* c32 */
25/* sin32 */
26/* cos32 */
27/* mpsin */
28/* mpcos */
29/* mpranred */
30/* mpsin1 */
31/* mpcos1 */
32/* */
33/* FILES NEEDED: endian.h mpa.h sincos32.h */
34/* mpa.c */
35/* */
36/* Multi Precision sin() and cos() function with p=32 for sin()*/
37/* cos() arcsin() and arccos() routines */
38/* In addition mpranred() routine performs range reduction of */
39/* a double number x into multi precision number y, */
40/* such that y=x-n*pi/2, abs(y)<pi/4, n=0,+-1,+-2,.... */
41/****************************************************************/
42#include "endian.h"
43#include "mpa.h"
44#include "sincos32.h"
45#include <math.h>
46#include <math_private.h>
47#include <stap-probe.h>
48
49#ifndef SECTION
50# define SECTION
51#endif
52
53/* Compute Multi-Precision sin() function for given p. Receive Multi Precision
54 number x and result stored at y. */
55static void
56SECTION
57ss32 (mp_no *x, mp_no *y, int p)
58{
59 int i;
60 double a;
61 mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
62 for (i = 1; i <= p; i++)
63 mpk.d[i] = 0;
64
65 __sqr (x, &x2, p);
66 __cpy (&oofac27, &gor, p);
67 __cpy (&gor, &sum, p);
68 for (a = 27.0; a > 1.0; a -= 2.0)
69 {
70 mpk.d[1] = a * (a - 1.0);
71 __mul (&gor, &mpk, &mpt1, p);
72 __cpy (&mpt1, &gor, p);
73 __mul (&x2, &sum, &mpt1, p);
74 __sub (&gor, &mpt1, &sum, p);
75 }
76 __mul (x, &sum, y, p);
77}
78
79/* Compute Multi-Precision cos() function for given p. Receive Multi Precision
80 number x and result stored at y. */
81static void
82SECTION
83cc32 (mp_no *x, mp_no *y, int p)
84{
85 int i;
86 double a;
87 mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}};
88 for (i = 1; i <= p; i++)
89 mpk.d[i] = 0;
90
91 __sqr (x, &x2, p);
92 mpk.d[1] = 27.0;
93 __mul (&oofac27, &mpk, &gor, p);
94 __cpy (&gor, &sum, p);
95 for (a = 26.0; a > 2.0; a -= 2.0)
96 {
97 mpk.d[1] = a * (a - 1.0);
98 __mul (&gor, &mpk, &mpt1, p);
99 __cpy (&mpt1, &gor, p);
100 __mul (&x2, &sum, &mpt1, p);
101 __sub (&gor, &mpt1, &sum, p);
102 }
103 __mul (&x2, &sum, y, p);
104}
105
106/* Compute both sin(x), cos(x) as Multi precision numbers. */
107void
108SECTION
109__c32 (mp_no *x, mp_no *y, mp_no *z, int p)
110{
111 mp_no u, t, t1, t2, c, s;
112 int i;
113 __cpy (x, &u, p);
114 u.e = u.e - 1;
115 cc32 (&u, &c, p);
116 ss32 (&u, &s, p);
117 for (i = 0; i < 24; i++)
118 {
119 __mul (&c, &s, &t, p);
120 __sub (&s, &t, &t1, p);
121 __add (&t1, &t1, &s, p);
122 __sub (&__mptwo, &c, &t1, p);
123 __mul (&t1, &c, &t2, p);
124 __add (&t2, &t2, &c, p);
125 }
126 __sub (&__mpone, &c, y, p);
127 __cpy (&s, z, p);
128}
129
130/* Receive double x and two double results of sin(x) and return result which is
131 more accurate, computing sin(x) with multi precision routine c32. */
132double
133SECTION
134__sin32 (double x, double res, double res1)
135{
136 int p;
137 mp_no a, b, c;
138 p = 32;
139 __dbl_mp (res, &a, p);
140 __dbl_mp (0.5 * (res1 - res), &b, p);
141 __add (&a, &b, &c, p);
142 if (x > 0.8)
143 {
144 __sub (&hp, &c, &a, p);
145 __c32 (&a, &b, &c, p);
146 }
147 else
148 __c32 (&c, &a, &b, p); /* b=sin(0.5*(res+res1)) */
149 __dbl_mp (x, &c, p); /* c = x */
150 __sub (&b, &c, &a, p);
151 /* if a > 0 return min (res, res1), otherwise return max (res, res1). */
152 if ((a.d[0] > 0 && res >= res1) || (a.d[0] <= 0 && res <= res1))
153 res = res1;
154 LIBC_PROBE (slowasin, 2, &res, &x);
155 return res;
156}
157
158/* Receive double x and two double results of cos(x) and return result which is
159 more accurate, computing cos(x) with multi precision routine c32. */
160double
161SECTION
162__cos32 (double x, double res, double res1)
163{
164 int p;
165 mp_no a, b, c;
166 p = 32;
167 __dbl_mp (res, &a, p);
168 __dbl_mp (0.5 * (res1 - res), &b, p);
169 __add (&a, &b, &c, p);
170 if (x > 2.4)
171 {
172 __sub (&pi, &c, &a, p);
173 __c32 (&a, &b, &c, p);
174 b.d[0] = -b.d[0];
175 }
176 else if (x > 0.8)
177 {
178 __sub (&hp, &c, &a, p);
179 __c32 (&a, &c, &b, p);
180 }
181 else
182 __c32 (&c, &b, &a, p); /* b=cos(0.5*(res+res1)) */
183 __dbl_mp (x, &c, p); /* c = x */
184 __sub (&b, &c, &a, p);
185 /* if a > 0 return max (res, res1), otherwise return min (res, res1). */
186 if ((a.d[0] > 0 && res <= res1) || (a.d[0] <= 0 && res >= res1))
187 res = res1;
188 LIBC_PROBE (slowacos, 2, &res, &x);
189 return res;
190}
191
192/* Compute sin() of double-length number (X + DX) as Multi Precision number and
193 return result as double. If REDUCE_RANGE is true, X is assumed to be the
194 original input and DX is ignored. */
195double
196SECTION
197__mpsin (double x, double dx, bool reduce_range)
198{
199 double y;
200 mp_no a, b, c, s;
201 int n;
202 int p = 32;
203
204 if (reduce_range)
205 {
206 n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */
207 __c32 (&a, &c, &s, p);
208 }
209 else
210 {
211 n = -1;
212 __dbl_mp (x, &b, p);
213 __dbl_mp (dx, &c, p);
214 __add (&b, &c, &a, p);
215 if (x > 0.8)
216 {
217 __sub (&hp, &a, &b, p);
218 __c32 (&b, &s, &c, p);
219 }
220 else
221 __c32 (&a, &c, &s, p); /* b = sin(x+dx) */
222 }
223
224 /* Convert result based on which quarter of unit circle y is in. */
225 switch (n)
226 {
227 case 1:
228 __mp_dbl (&c, &y, p);
229 break;
230
231 case 3:
232 __mp_dbl (&c, &y, p);
233 y = -y;
234 break;
235
236 case 2:
237 __mp_dbl (&s, &y, p);
238 y = -y;
239 break;
240
241 /* Quadrant not set, so the result must be sin (X + DX), which is also in
242 S. */
243 case 0:
244 default:
245 __mp_dbl (&s, &y, p);
246 }
247 LIBC_PROBE (slowsin, 3, &x, &dx, &y);
248 return y;
249}
250
251/* Compute cos() of double-length number (X + DX) as Multi Precision number and
252 return result as double. If REDUCE_RANGE is true, X is assumed to be the
253 original input and DX is ignored. */
254double
255SECTION
256__mpcos (double x, double dx, bool reduce_range)
257{
258 double y;
259 mp_no a, b, c, s;
260 int n;
261 int p = 32;
262
263 if (reduce_range)
264 {
265 n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */
266 __c32 (&a, &c, &s, p);
267 }
268 else
269 {
270 n = -1;
271 __dbl_mp (x, &b, p);
272 __dbl_mp (dx, &c, p);
273 __add (&b, &c, &a, p);
274 if (x > 0.8)
275 {
276 __sub (&hp, &a, &b, p);
277 __c32 (&b, &s, &c, p);
278 }
279 else
280 __c32 (&a, &c, &s, p); /* a = cos(x+dx) */
281 }
282
283 /* Convert result based on which quarter of unit circle y is in. */
284 switch (n)
285 {
286 case 1:
287 __mp_dbl (&s, &y, p);
288 y = -y;
289 break;
290
291 case 3:
292 __mp_dbl (&s, &y, p);
293 break;
294
295 case 2:
296 __mp_dbl (&c, &y, p);
297 y = -y;
298 break;
299
300 /* Quadrant not set, so the result must be cos (X + DX), which is also
301 stored in C. */
302 case 0:
303 default:
304 __mp_dbl (&c, &y, p);
305 }
306 LIBC_PROBE (slowcos, 3, &x, &dx, &y);
307 return y;
308}
309
310/* Perform range reduction of a double number x into multi precision number y,
311 such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ...
312 Return int which indicates in which quarter of circle x is. */
313int
314SECTION
315__mpranred (double x, mp_no *y, int p)
316{
317 number v;
318 double t, xn;
319 int i, k, n;
320 mp_no a, b, c;
321
322 if (fabs (x) < 2.8e14)
323 {
324 t = (x * hpinv.d + toint.d);
325 xn = t - toint.d;
326 v.d = t;
327 n = v.i[LOW_HALF] & 3;
328 __dbl_mp (xn, &a, p);
329 __mul (&a, &hp, &b, p);
330 __dbl_mp (x, &c, p);
331 __sub (&c, &b, y, p);
332 return n;
333 }
334 else
335 {
336 /* If x is very big more precision required. */
337 __dbl_mp (x, &a, p);
338 a.d[0] = 1.0;
339 k = a.e - 5;
340 if (k < 0)
341 k = 0;
342 b.e = -k;
343 b.d[0] = 1.0;
344 for (i = 0; i < p; i++)
345 b.d[i + 1] = toverp[i + k];
346 __mul (&a, &b, &c, p);
347 t = c.d[c.e];
348 for (i = 1; i <= p - c.e; i++)
349 c.d[i] = c.d[i + c.e];
350 for (i = p + 1 - c.e; i <= p; i++)
351 c.d[i] = 0;
352 c.e = 0;
353 if (c.d[1] >= HALFRAD)
354 {
355 t += 1.0;
356 __sub (&c, &__mpone, &b, p);
357 __mul (&b, &hp, y, p);
358 }
359 else
360 __mul (&c, &hp, y, p);
361 n = (int) t;
362 if (x < 0)
363 {
364 y->d[0] = -y->d[0];
365 n = -n;
366 }
367 return (n & 3);
368 }
369}
370