]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 | 3 | * written by International Business Machines Corp. |
f7a9f785 | 4 | * Copyright (C) 2001-2016 Free Software Foundation, Inc. |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
50944bca | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
e4d82761 UD |
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" | |
0e9be4db | 45 | #include <math.h> |
1ed0291c | 46 | #include <math_private.h> |
f3fd2628 | 47 | #include <stap-probe.h> |
e4d82761 | 48 | |
31d3cc00 UD |
49 | #ifndef SECTION |
50 | # define SECTION | |
51 | #endif | |
52 | ||
97a0650b SP |
53 | /* Compute Multi-Precision sin() function for given p. Receive Multi Precision |
54 | number x and result stored at y. */ | |
31d3cc00 UD |
55 | static void |
56 | SECTION | |
97a0650b SP |
57 | ss32 (mp_no *x, mp_no *y, int p) |
58 | { | |
e4d82761 | 59 | int i; |
50944bca | 60 | double a; |
97a0650b SP |
61 | mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; |
62 | for (i = 1; i <= p; i++) | |
63 | mpk.d[i] = 0; | |
e4d82761 | 64 | |
97a0650b SP |
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); | |
e4d82761 UD |
77 | } |
78 | ||
97a0650b SP |
79 | /* Compute Multi-Precision cos() function for given p. Receive Multi Precision |
80 | number x and result stored at y. */ | |
31d3cc00 UD |
81 | static void |
82 | SECTION | |
97a0650b SP |
83 | cc32 (mp_no *x, mp_no *y, int p) |
84 | { | |
e4d82761 | 85 | int i; |
50944bca | 86 | double a; |
97a0650b SP |
87 | mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; |
88 | for (i = 1; i <= p; i++) | |
89 | mpk.d[i] = 0; | |
e4d82761 | 90 | |
97a0650b SP |
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); | |
e4d82761 UD |
104 | } |
105 | ||
97a0650b | 106 | /* Compute both sin(x), cos(x) as Multi precision numbers. */ |
31d3cc00 UD |
107 | void |
108 | SECTION | |
97a0650b SP |
109 | __c32 (mp_no *x, mp_no *y, mp_no *z, int p) |
110 | { | |
111 | mp_no u, t, t1, t2, c, s; | |
e4d82761 | 112 | int i; |
97a0650b SP |
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); | |
107a5bf0 | 122 | __sub (&__mptwo, &c, &t1, p); |
97a0650b SP |
123 | __mul (&t1, &c, &t2, p); |
124 | __add (&t2, &t2, &c, p); | |
125 | } | |
107a5bf0 | 126 | __sub (&__mpone, &c, y, p); |
97a0650b | 127 | __cpy (&s, z, p); |
e4d82761 UD |
128 | } |
129 | ||
97a0650b SP |
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. */ | |
31d3cc00 UD |
132 | double |
133 | SECTION | |
97a0650b SP |
134 | __sin32 (double x, double res, double res1) |
135 | { | |
e4d82761 | 136 | int p; |
97a0650b SP |
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). */ | |
c79a1204 SP |
152 | if ((a.d[0] > 0 && res >= res1) || (a.d[0] <= 0 && res <= res1)) |
153 | res = res1; | |
f3fd2628 | 154 | LIBC_PROBE (slowasin, 2, &res, &x); |
c79a1204 | 155 | return res; |
e4d82761 UD |
156 | } |
157 | ||
97a0650b SP |
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. */ | |
31d3cc00 UD |
160 | double |
161 | SECTION | |
97a0650b SP |
162 | __cos32 (double x, double res, double res1) |
163 | { | |
e4d82761 | 164 | int p; |
97a0650b SP |
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). */ | |
c79a1204 SP |
186 | if ((a.d[0] > 0 && res <= res1) || (a.d[0] <= 0 && res >= res1)) |
187 | res = res1; | |
f3fd2628 | 188 | LIBC_PROBE (slowacos, 2, &res, &x); |
c79a1204 | 189 | return res; |
e4d82761 UD |
190 | } |
191 | ||
09544cbc SP |
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. */ | |
31d3cc00 UD |
195 | double |
196 | SECTION | |
09544cbc | 197 | __mpsin (double x, double dx, bool reduce_range) |
97a0650b | 198 | { |
e4d82761 | 199 | double y; |
09544cbc SP |
200 | mp_no a, b, c, s; |
201 | int n; | |
202 | int p = 32; | |
203 | ||
204 | if (reduce_range) | |
97a0650b | 205 | { |
09544cbc SP |
206 | n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ |
207 | __c32 (&a, &c, &s, p); | |
97a0650b SP |
208 | } |
209 | else | |
09544cbc SP |
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 | } | |
f3fd2628 | 247 | LIBC_PROBE (slowsin, 3, &x, &dx, &y); |
e4d82761 UD |
248 | return y; |
249 | } | |
250 | ||
09544cbc SP |
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. */ | |
31d3cc00 UD |
254 | double |
255 | SECTION | |
09544cbc | 256 | __mpcos (double x, double dx, bool reduce_range) |
97a0650b | 257 | { |
e4d82761 | 258 | double y; |
09544cbc SP |
259 | mp_no a, b, c, s; |
260 | int n; | |
261 | int p = 32; | |
262 | ||
263 | if (reduce_range) | |
97a0650b | 264 | { |
09544cbc SP |
265 | n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ |
266 | __c32 (&a, &c, &s, p); | |
97a0650b SP |
267 | } |
268 | else | |
09544cbc SP |
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 | } | |
f3fd2628 | 306 | LIBC_PROBE (slowcos, 3, &x, &dx, &y); |
e4d82761 UD |
307 | return y; |
308 | } | |
309 | ||
97a0650b SP |
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. */ | |
31d3cc00 UD |
313 | int |
314 | SECTION | |
97a0650b | 315 | __mpranred (double x, mp_no *y, int p) |
e4d82761 UD |
316 | { |
317 | number v; | |
97a0650b SP |
318 | double t, xn; |
319 | int i, k, n; | |
320 | mp_no a, b, c; | |
50944bca | 321 | |
0e9be4db | 322 | if (fabs (x) < 2.8e14) |
97a0650b SP |
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; | |
107a5bf0 | 356 | __sub (&c, &__mpone, &b, p); |
97a0650b SP |
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); | |
e4d82761 | 368 | } |
e4d82761 | 369 | } |