2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
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.
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.
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/>.
19 /****************************************************************/
20 /* MODULE_NAME: sincos32.c */
33 /* FILES NEEDED: endian.h mpa.h sincos32.h */
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 /****************************************************************/
46 #include <math_private.h>
47 #include <stap-probe.h>
53 /* Compute Multi-Precision sin() function for given p. Receive Multi Precision
54 number x and result stored at y. */
57 ss32 (mp_no
*x
, mp_no
*y
, int p
)
61 mp_no mpt1
, x2
, gor
, sum
, mpk
= {1, {1.0}};
62 for (i
= 1; i
<= p
; i
++)
66 __cpy (&oofac27
, &gor
, p
);
67 __cpy (&gor
, &sum
, p
);
68 for (a
= 27.0; a
> 1.0; a
-= 2.0)
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
);
76 __mul (x
, &sum
, y
, p
);
79 /* Compute Multi-Precision cos() function for given p. Receive Multi Precision
80 number x and result stored at y. */
83 cc32 (mp_no
*x
, mp_no
*y
, int p
)
87 mp_no mpt1
, x2
, gor
, sum
, mpk
= {1, {1.0}};
88 for (i
= 1; i
<= p
; i
++)
93 __mul (&oofac27
, &mpk
, &gor
, p
);
94 __cpy (&gor
, &sum
, p
);
95 for (a
= 26.0; a
> 2.0; a
-= 2.0)
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
);
103 __mul (&x2
, &sum
, y
, p
);
106 /* Compute both sin(x), cos(x) as Multi precision numbers. */
109 __c32 (mp_no
*x
, mp_no
*y
, mp_no
*z
, int p
)
111 mp_no u
, t
, t1
, t2
, c
, s
;
117 for (i
= 0; i
< 24; i
++)
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
);
126 __sub (&__mpone
, &c
, y
, p
);
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. */
134 __sin32 (double x
, double res
, double res1
)
139 __dbl_mp (res
, &a
, p
);
140 __dbl_mp (0.5 * (res1
- res
), &b
, p
);
141 __add (&a
, &b
, &c
, p
);
144 __sub (&hp
, &c
, &a
, p
);
145 __c32 (&a
, &b
, &c
, p
);
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
))
154 LIBC_PROBE (slowasin
, 2, &res
, &x
);
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. */
162 __cos32 (double x
, double res
, double res1
)
167 __dbl_mp (res
, &a
, p
);
168 __dbl_mp (0.5 * (res1
- res
), &b
, p
);
169 __add (&a
, &b
, &c
, p
);
172 __sub (&pi
, &c
, &a
, p
);
173 __c32 (&a
, &b
, &c
, p
);
178 __sub (&hp
, &c
, &a
, p
);
179 __c32 (&a
, &c
, &b
, p
);
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
))
188 LIBC_PROBE (slowacos
, 2, &res
, &x
);
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. */
197 __mpsin (double x
, double dx
, bool reduce_range
)
206 n
= __mpranred (x
, &a
, p
); /* n is 0, 1, 2 or 3. */
207 __c32 (&a
, &c
, &s
, p
);
213 __dbl_mp (dx
, &c
, p
);
214 __add (&b
, &c
, &a
, p
);
217 __sub (&hp
, &a
, &b
, p
);
218 __c32 (&b
, &s
, &c
, p
);
221 __c32 (&a
, &c
, &s
, p
); /* b = sin(x+dx) */
224 /* Convert result based on which quarter of unit circle y is in. */
228 __mp_dbl (&c
, &y
, p
);
232 __mp_dbl (&c
, &y
, p
);
237 __mp_dbl (&s
, &y
, p
);
241 /* Quadrant not set, so the result must be sin (X + DX), which is also in
245 __mp_dbl (&s
, &y
, p
);
247 LIBC_PROBE (slowsin
, 3, &x
, &dx
, &y
);
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. */
256 __mpcos (double x
, double dx
, bool reduce_range
)
265 n
= __mpranred (x
, &a
, p
); /* n is 0, 1, 2 or 3. */
266 __c32 (&a
, &c
, &s
, p
);
272 __dbl_mp (dx
, &c
, p
);
273 __add (&b
, &c
, &a
, p
);
276 __sub (&hp
, &a
, &b
, p
);
277 __c32 (&b
, &s
, &c
, p
);
280 __c32 (&a
, &c
, &s
, p
); /* a = cos(x+dx) */
283 /* Convert result based on which quarter of unit circle y is in. */
287 __mp_dbl (&s
, &y
, p
);
292 __mp_dbl (&s
, &y
, p
);
296 __mp_dbl (&c
, &y
, p
);
300 /* Quadrant not set, so the result must be cos (X + DX), which is also
304 __mp_dbl (&c
, &y
, p
);
306 LIBC_PROBE (slowcos
, 3, &x
, &dx
, &y
);
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. */
315 __mpranred (double x
, mp_no
*y
, int p
)
322 if (fabs (x
) < 2.8e14
)
324 t
= (x
* hpinv
.d
+ toint
.d
);
327 n
= v
.i
[LOW_HALF
] & 3;
328 __dbl_mp (xn
, &a
, p
);
329 __mul (&a
, &hp
, &b
, p
);
331 __sub (&c
, &b
, y
, p
);
336 /* If x is very big more precision required. */
344 for (i
= 0; i
< p
; i
++)
345 b
.d
[i
+ 1] = toverp
[i
+ k
];
346 __mul (&a
, &b
, &c
, p
);
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
++)
353 if (c
.d
[1] >= HALFRAD
)
356 __sub (&c
, &__mpone
, &b
, p
);
357 __mul (&b
, &hp
, y
, p
);
360 __mul (&c
, &hp
, y
, p
);