1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 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 <http://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. */ |
55 | static void |
56 | SECTION |
57 | ss32 (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. */ |
81 | static void |
82 | SECTION |
83 | cc32 (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. */ |
107 | void |
108 | SECTION |
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. */ |
132 | double |
133 | SECTION |
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. */ |
160 | double |
161 | SECTION |
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. */ |
195 | double |
196 | SECTION |
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. */ |
254 | double |
255 | SECTION |
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. */ |
313 | int |
314 | SECTION |
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 | |