1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2021 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/* Compute sin() of double-length number (X + DX) as Multi Precision number and
131 return result as double. If REDUCE_RANGE is true, X is assumed to be the
132 original input and DX is ignored. */
133double
134SECTION
135__mpsin (double x, double dx, bool reduce_range)
136{
137 double y;
138 mp_no a, b, c, s;
139 int n;
140 int p = 32;
141
142 if (reduce_range)
143 {
144 n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */
145 __c32 (&a, &c, &s, p);
146 }
147 else
148 {
149 n = -1;
150 __dbl_mp (x, &b, p);
151 __dbl_mp (dx, &c, p);
152 __add (&b, &c, &a, p);
153 if (x > 0.8)
154 {
155 __sub (&hp, &a, &b, p);
156 __c32 (&b, &s, &c, p);
157 }
158 else
159 __c32 (&a, &c, &s, p); /* b = sin(x+dx) */
160 }
161
162 /* Convert result based on which quarter of unit circle y is in. */
163 switch (n)
164 {
165 case 1:
166 __mp_dbl (&c, &y, p);
167 break;
168
169 case 3:
170 __mp_dbl (&c, &y, p);
171 y = -y;
172 break;
173
174 case 2:
175 __mp_dbl (&s, &y, p);
176 y = -y;
177 break;
178
179 /* Quadrant not set, so the result must be sin (X + DX), which is also in
180 S. */
181 case 0:
182 default:
183 __mp_dbl (&s, &y, p);
184 }
185 LIBC_PROBE (slowsin, 3, &x, &dx, &y);
186 return y;
187}
188
189/* Compute cos() of double-length number (X + DX) as Multi Precision number and
190 return result as double. If REDUCE_RANGE is true, X is assumed to be the
191 original input and DX is ignored. */
192double
193SECTION
194__mpcos (double x, double dx, bool reduce_range)
195{
196 double y;
197 mp_no a, b, c, s;
198 int n;
199 int p = 32;
200
201 if (reduce_range)
202 {
203 n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */
204 __c32 (&a, &c, &s, p);
205 }
206 else
207 {
208 n = -1;
209 __dbl_mp (x, &b, p);
210 __dbl_mp (dx, &c, p);
211 __add (&b, &c, &a, p);
212 if (x > 0.8)
213 {
214 __sub (&hp, &a, &b, p);
215 __c32 (&b, &s, &c, p);
216 }
217 else
218 __c32 (&a, &c, &s, p); /* a = cos(x+dx) */
219 }
220
221 /* Convert result based on which quarter of unit circle y is in. */
222 switch (n)
223 {
224 case 1:
225 __mp_dbl (&s, &y, p);
226 y = -y;
227 break;
228
229 case 3:
230 __mp_dbl (&s, &y, p);
231 break;
232
233 case 2:
234 __mp_dbl (&c, &y, p);
235 y = -y;
236 break;
237
238 /* Quadrant not set, so the result must be cos (X + DX), which is also
239 stored in C. */
240 case 0:
241 default:
242 __mp_dbl (&c, &y, p);
243 }
244 LIBC_PROBE (slowcos, 3, &x, &dx, &y);
245 return y;
246}
247
248/* Perform range reduction of a double number x into multi precision number y,
249 such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ...
250 Return int which indicates in which quarter of circle x is. */
251int
252SECTION
253__mpranred (double x, mp_no *y, int p)
254{
255 number v;
256 double t, xn;
257 int i, k, n;
258 mp_no a, b, c;
259
260 if (fabs (x) < 2.8e14)
261 {
262 t = (x * hpinv.d + toint.d);
263 xn = t - toint.d;
264 v.d = t;
265 n = v.i[LOW_HALF] & 3;
266 __dbl_mp (xn, &a, p);
267 __mul (&a, &hp, &b, p);
268 __dbl_mp (x, &c, p);
269 __sub (&c, &b, y, p);
270 return n;
271 }
272 else
273 {
274 /* If x is very big more precision required. */
275 __dbl_mp (x, &a, p);
276 a.d[0] = 1.0;
277 k = a.e - 5;
278 if (k < 0)
279 k = 0;
280 b.e = -k;
281 b.d[0] = 1.0;
282 for (i = 0; i < p; i++)
283 b.d[i + 1] = toverp[i + k];
284 __mul (&a, &b, &c, p);
285 t = c.d[c.e];
286 for (i = 1; i <= p - c.e; i++)
287 c.d[i] = c.d[i + c.e];
288 for (i = p + 1 - c.e; i <= p; i++)
289 c.d[i] = 0;
290 c.e = 0;
291 if (c.d[1] >= HALFRAD)
292 {
293 t += 1.0;
294 __sub (&c, &__mpone, &b, p);
295 __mul (&b, &hp, y, p);
296 }
297 else
298 __mul (&c, &hp, y, p);
299 n = (int) t;
300 if (x < 0)
301 {
302 y->d[0] = -y->d[0];
303 n = -n;
304 }
305 return (n & 3);
306 }
307}
308