]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 | 3 | * written by International Business Machines Corp. |
b168057a | 4 | * Copyright (C) 2001-2015 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" | |
1ed0291c | 45 | #include <math_private.h> |
f3fd2628 | 46 | #include <stap-probe.h> |
e4d82761 | 47 | |
31d3cc00 UD |
48 | #ifndef SECTION |
49 | # define SECTION | |
50 | #endif | |
51 | ||
97a0650b SP |
52 | /* Compute Multi-Precision sin() function for given p. Receive Multi Precision |
53 | number x and result stored at y. */ | |
31d3cc00 UD |
54 | static void |
55 | SECTION | |
97a0650b SP |
56 | ss32 (mp_no *x, mp_no *y, int p) |
57 | { | |
e4d82761 | 58 | int i; |
50944bca | 59 | double a; |
97a0650b SP |
60 | mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; |
61 | for (i = 1; i <= p; i++) | |
62 | mpk.d[i] = 0; | |
e4d82761 | 63 | |
97a0650b SP |
64 | __sqr (x, &x2, p); |
65 | __cpy (&oofac27, &gor, p); | |
66 | __cpy (&gor, &sum, p); | |
67 | for (a = 27.0; a > 1.0; a -= 2.0) | |
68 | { | |
69 | mpk.d[1] = a * (a - 1.0); | |
70 | __mul (&gor, &mpk, &mpt1, p); | |
71 | __cpy (&mpt1, &gor, p); | |
72 | __mul (&x2, &sum, &mpt1, p); | |
73 | __sub (&gor, &mpt1, &sum, p); | |
74 | } | |
75 | __mul (x, &sum, y, p); | |
e4d82761 UD |
76 | } |
77 | ||
97a0650b SP |
78 | /* Compute Multi-Precision cos() function for given p. Receive Multi Precision |
79 | number x and result stored at y. */ | |
31d3cc00 UD |
80 | static void |
81 | SECTION | |
97a0650b SP |
82 | cc32 (mp_no *x, mp_no *y, int p) |
83 | { | |
e4d82761 | 84 | int i; |
50944bca | 85 | double a; |
97a0650b SP |
86 | mp_no mpt1, x2, gor, sum, mpk = {1, {1.0}}; |
87 | for (i = 1; i <= p; i++) | |
88 | mpk.d[i] = 0; | |
e4d82761 | 89 | |
97a0650b SP |
90 | __sqr (x, &x2, p); |
91 | mpk.d[1] = 27.0; | |
92 | __mul (&oofac27, &mpk, &gor, p); | |
93 | __cpy (&gor, &sum, p); | |
94 | for (a = 26.0; a > 2.0; a -= 2.0) | |
95 | { | |
96 | mpk.d[1] = a * (a - 1.0); | |
97 | __mul (&gor, &mpk, &mpt1, p); | |
98 | __cpy (&mpt1, &gor, p); | |
99 | __mul (&x2, &sum, &mpt1, p); | |
100 | __sub (&gor, &mpt1, &sum, p); | |
101 | } | |
102 | __mul (&x2, &sum, y, p); | |
e4d82761 UD |
103 | } |
104 | ||
97a0650b | 105 | /* Compute both sin(x), cos(x) as Multi precision numbers. */ |
31d3cc00 UD |
106 | void |
107 | SECTION | |
97a0650b SP |
108 | __c32 (mp_no *x, mp_no *y, mp_no *z, int p) |
109 | { | |
110 | mp_no u, t, t1, t2, c, s; | |
e4d82761 | 111 | int i; |
97a0650b SP |
112 | __cpy (x, &u, p); |
113 | u.e = u.e - 1; | |
114 | cc32 (&u, &c, p); | |
115 | ss32 (&u, &s, p); | |
116 | for (i = 0; i < 24; i++) | |
117 | { | |
118 | __mul (&c, &s, &t, p); | |
119 | __sub (&s, &t, &t1, p); | |
120 | __add (&t1, &t1, &s, p); | |
107a5bf0 | 121 | __sub (&__mptwo, &c, &t1, p); |
97a0650b SP |
122 | __mul (&t1, &c, &t2, p); |
123 | __add (&t2, &t2, &c, p); | |
124 | } | |
107a5bf0 | 125 | __sub (&__mpone, &c, y, p); |
97a0650b | 126 | __cpy (&s, z, p); |
e4d82761 UD |
127 | } |
128 | ||
97a0650b SP |
129 | /* Receive double x and two double results of sin(x) and return result which is |
130 | more accurate, computing sin(x) with multi precision routine c32. */ | |
31d3cc00 UD |
131 | double |
132 | SECTION | |
97a0650b SP |
133 | __sin32 (double x, double res, double res1) |
134 | { | |
e4d82761 | 135 | int p; |
97a0650b SP |
136 | mp_no a, b, c; |
137 | p = 32; | |
138 | __dbl_mp (res, &a, p); | |
139 | __dbl_mp (0.5 * (res1 - res), &b, p); | |
140 | __add (&a, &b, &c, p); | |
141 | if (x > 0.8) | |
142 | { | |
143 | __sub (&hp, &c, &a, p); | |
144 | __c32 (&a, &b, &c, p); | |
145 | } | |
146 | else | |
147 | __c32 (&c, &a, &b, p); /* b=sin(0.5*(res+res1)) */ | |
148 | __dbl_mp (x, &c, p); /* c = x */ | |
149 | __sub (&b, &c, &a, p); | |
150 | /* if a > 0 return min (res, res1), otherwise return max (res, res1). */ | |
c79a1204 SP |
151 | if ((a.d[0] > 0 && res >= res1) || (a.d[0] <= 0 && res <= res1)) |
152 | res = res1; | |
f3fd2628 | 153 | LIBC_PROBE (slowasin, 2, &res, &x); |
c79a1204 | 154 | return res; |
e4d82761 UD |
155 | } |
156 | ||
97a0650b SP |
157 | /* Receive double x and two double results of cos(x) and return result which is |
158 | more accurate, computing cos(x) with multi precision routine c32. */ | |
31d3cc00 UD |
159 | double |
160 | SECTION | |
97a0650b SP |
161 | __cos32 (double x, double res, double res1) |
162 | { | |
e4d82761 | 163 | int p; |
97a0650b SP |
164 | mp_no a, b, c; |
165 | p = 32; | |
166 | __dbl_mp (res, &a, p); | |
167 | __dbl_mp (0.5 * (res1 - res), &b, p); | |
168 | __add (&a, &b, &c, p); | |
169 | if (x > 2.4) | |
170 | { | |
171 | __sub (&pi, &c, &a, p); | |
172 | __c32 (&a, &b, &c, p); | |
173 | b.d[0] = -b.d[0]; | |
174 | } | |
175 | else if (x > 0.8) | |
176 | { | |
177 | __sub (&hp, &c, &a, p); | |
178 | __c32 (&a, &c, &b, p); | |
179 | } | |
180 | else | |
181 | __c32 (&c, &b, &a, p); /* b=cos(0.5*(res+res1)) */ | |
182 | __dbl_mp (x, &c, p); /* c = x */ | |
183 | __sub (&b, &c, &a, p); | |
184 | /* if a > 0 return max (res, res1), otherwise return min (res, res1). */ | |
c79a1204 SP |
185 | if ((a.d[0] > 0 && res <= res1) || (a.d[0] <= 0 && res >= res1)) |
186 | res = res1; | |
f3fd2628 | 187 | LIBC_PROBE (slowacos, 2, &res, &x); |
c79a1204 | 188 | return res; |
e4d82761 UD |
189 | } |
190 | ||
09544cbc SP |
191 | /* Compute sin() of double-length number (X + DX) as Multi Precision number and |
192 | return result as double. If REDUCE_RANGE is true, X is assumed to be the | |
193 | original input and DX is ignored. */ | |
31d3cc00 UD |
194 | double |
195 | SECTION | |
09544cbc | 196 | __mpsin (double x, double dx, bool reduce_range) |
97a0650b | 197 | { |
e4d82761 | 198 | double y; |
09544cbc SP |
199 | mp_no a, b, c, s; |
200 | int n; | |
201 | int p = 32; | |
202 | ||
203 | if (reduce_range) | |
97a0650b | 204 | { |
09544cbc SP |
205 | n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ |
206 | __c32 (&a, &c, &s, p); | |
97a0650b SP |
207 | } |
208 | else | |
09544cbc SP |
209 | { |
210 | n = -1; | |
211 | __dbl_mp (x, &b, p); | |
212 | __dbl_mp (dx, &c, p); | |
213 | __add (&b, &c, &a, p); | |
214 | if (x > 0.8) | |
215 | { | |
216 | __sub (&hp, &a, &b, p); | |
217 | __c32 (&b, &s, &c, p); | |
218 | } | |
219 | else | |
220 | __c32 (&a, &c, &s, p); /* b = sin(x+dx) */ | |
221 | } | |
222 | ||
223 | /* Convert result based on which quarter of unit circle y is in. */ | |
224 | switch (n) | |
225 | { | |
226 | case 1: | |
227 | __mp_dbl (&c, &y, p); | |
228 | break; | |
229 | ||
230 | case 3: | |
231 | __mp_dbl (&c, &y, p); | |
232 | y = -y; | |
233 | break; | |
234 | ||
235 | case 2: | |
236 | __mp_dbl (&s, &y, p); | |
237 | y = -y; | |
238 | break; | |
239 | ||
240 | /* Quadrant not set, so the result must be sin (X + DX), which is also in | |
241 | S. */ | |
242 | case 0: | |
243 | default: | |
244 | __mp_dbl (&s, &y, p); | |
245 | } | |
f3fd2628 | 246 | LIBC_PROBE (slowsin, 3, &x, &dx, &y); |
e4d82761 UD |
247 | return y; |
248 | } | |
249 | ||
09544cbc SP |
250 | /* Compute cos() of double-length number (X + DX) as Multi Precision number and |
251 | return result as double. If REDUCE_RANGE is true, X is assumed to be the | |
252 | original input and DX is ignored. */ | |
31d3cc00 UD |
253 | double |
254 | SECTION | |
09544cbc | 255 | __mpcos (double x, double dx, bool reduce_range) |
97a0650b | 256 | { |
e4d82761 | 257 | double y; |
09544cbc SP |
258 | mp_no a, b, c, s; |
259 | int n; | |
260 | int p = 32; | |
261 | ||
262 | if (reduce_range) | |
97a0650b | 263 | { |
09544cbc SP |
264 | n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ |
265 | __c32 (&a, &c, &s, p); | |
97a0650b SP |
266 | } |
267 | else | |
09544cbc SP |
268 | { |
269 | n = -1; | |
270 | __dbl_mp (x, &b, p); | |
271 | __dbl_mp (dx, &c, p); | |
272 | __add (&b, &c, &a, p); | |
273 | if (x > 0.8) | |
274 | { | |
275 | __sub (&hp, &a, &b, p); | |
276 | __c32 (&b, &s, &c, p); | |
277 | } | |
278 | else | |
279 | __c32 (&a, &c, &s, p); /* a = cos(x+dx) */ | |
280 | } | |
281 | ||
282 | /* Convert result based on which quarter of unit circle y is in. */ | |
283 | switch (n) | |
284 | { | |
285 | case 1: | |
286 | __mp_dbl (&s, &y, p); | |
287 | y = -y; | |
288 | break; | |
289 | ||
290 | case 3: | |
291 | __mp_dbl (&s, &y, p); | |
292 | break; | |
293 | ||
294 | case 2: | |
295 | __mp_dbl (&c, &y, p); | |
296 | y = -y; | |
297 | break; | |
298 | ||
299 | /* Quadrant not set, so the result must be cos (X + DX), which is also | |
300 | stored in C. */ | |
301 | case 0: | |
302 | default: | |
303 | __mp_dbl (&c, &y, p); | |
304 | } | |
f3fd2628 | 305 | LIBC_PROBE (slowcos, 3, &x, &dx, &y); |
e4d82761 UD |
306 | return y; |
307 | } | |
308 | ||
97a0650b SP |
309 | /* Perform range reduction of a double number x into multi precision number y, |
310 | such that y = x - n * pi / 2, abs (y) < pi / 4, n = 0, +-1, +-2, ... | |
311 | Return int which indicates in which quarter of circle x is. */ | |
31d3cc00 UD |
312 | int |
313 | SECTION | |
97a0650b | 314 | __mpranred (double x, mp_no *y, int p) |
e4d82761 UD |
315 | { |
316 | number v; | |
97a0650b SP |
317 | double t, xn; |
318 | int i, k, n; | |
319 | mp_no a, b, c; | |
50944bca | 320 | |
97a0650b SP |
321 | if (ABS (x) < 2.8e14) |
322 | { | |
323 | t = (x * hpinv.d + toint.d); | |
324 | xn = t - toint.d; | |
325 | v.d = t; | |
326 | n = v.i[LOW_HALF] & 3; | |
327 | __dbl_mp (xn, &a, p); | |
328 | __mul (&a, &hp, &b, p); | |
329 | __dbl_mp (x, &c, p); | |
330 | __sub (&c, &b, y, p); | |
331 | return n; | |
332 | } | |
333 | else | |
334 | { | |
335 | /* If x is very big more precision required. */ | |
336 | __dbl_mp (x, &a, p); | |
337 | a.d[0] = 1.0; | |
338 | k = a.e - 5; | |
339 | if (k < 0) | |
340 | k = 0; | |
341 | b.e = -k; | |
342 | b.d[0] = 1.0; | |
343 | for (i = 0; i < p; i++) | |
344 | b.d[i + 1] = toverp[i + k]; | |
345 | __mul (&a, &b, &c, p); | |
346 | t = c.d[c.e]; | |
347 | for (i = 1; i <= p - c.e; i++) | |
348 | c.d[i] = c.d[i + c.e]; | |
349 | for (i = p + 1 - c.e; i <= p; i++) | |
350 | c.d[i] = 0; | |
351 | c.e = 0; | |
352 | if (c.d[1] >= HALFRAD) | |
353 | { | |
354 | t += 1.0; | |
107a5bf0 | 355 | __sub (&c, &__mpone, &b, p); |
97a0650b SP |
356 | __mul (&b, &hp, y, p); |
357 | } | |
358 | else | |
359 | __mul (&c, &hp, y, p); | |
360 | n = (int) t; | |
361 | if (x < 0) | |
362 | { | |
363 | y->d[0] = -y->d[0]; | |
364 | n = -n; | |
365 | } | |
366 | return (n & 3); | |
e4d82761 | 367 | } |
e4d82761 | 368 | } |