2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001, 2011 Free Software Foundation
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 /****************************************************************/
45 #include "math_private.h"
51 /****************************************************************/
52 /* Compute Multi-Precision sin() function for given p. Receive */
53 /* Multi Precision number x and result stored at y */
54 /****************************************************************/
57 ss32(mp_no
*x
, mp_no
*y
, int p
) {
62 static const mp_no mpone
= {1,{1.0,1.0}};
64 mp_no mpt1
,x2
,gor
,sum
,mpk
={1,{1.0}};
68 for (i
=1;i
<=p
;i
++) mpk
.d
[i
]=0;
71 __cpy(&oofac27
,&gor
,p
);
73 for (a
=27.0;a
>1.0;a
-=2.0) {
75 __mul(&gor
,&mpk
,&mpt1
,p
);
77 __mul(&x2
,&sum
,&mpt1
,p
);
78 __sub(&gor
,&mpt1
,&sum
,p
);
83 /**********************************************************************/
84 /* Compute Multi-Precision cos() function for given p. Receive Multi */
85 /* Precision number x and result stored at y */
86 /**********************************************************************/
89 cc32(mp_no
*x
, mp_no
*y
, int p
) {
94 static const mp_no mpone
= {1,{1.0,1.0}};
96 mp_no mpt1
,x2
,gor
,sum
,mpk
={1,{1.0}};
100 for (i
=1;i
<=p
;i
++) mpk
.d
[i
]=0;
104 __mul(&oofac27
,&mpk
,&gor
,p
);
106 for (a
=26.0;a
>2.0;a
-=2.0) {
108 __mul(&gor
,&mpk
,&mpt1
,p
);
110 __mul(&x2
,&sum
,&mpt1
,p
);
111 __sub(&gor
,&mpt1
,&sum
,p
);
116 /***************************************************************************/
117 /* c32() computes both sin(x), cos(x) as Multi precision numbers */
118 /***************************************************************************/
121 __c32(mp_no
*x
, mp_no
*y
, mp_no
*z
, int p
) {
122 static const mp_no mpt
={1,{1.0,2.0}}, one
={1,{1.0,1.0}};
133 __sub(&mpt
,&c
,&t1
,p
);
141 /************************************************************************/
142 /*Routine receive double x and two double results of sin(x) and return */
143 /*result which is more accurate */
144 /*Computing sin(x) with multi precision routine c32 */
145 /************************************************************************/
148 __sin32(double x
, double res
, double res1
) {
153 __dbl_mp(0.5*(res1
-res
),&b
,p
);
156 { __sub(&hp
,&c
,&a
,p
);
159 else __c32(&c
,&a
,&b
,p
); /* b=sin(0.5*(res+res1)) */
160 __dbl_mp(x
,&c
,p
); /* c = x */
162 /* if a>0 return min(res,res1), otherwise return max(res,res1) */
163 if (a
.d
[0]>0) return (res
<res1
)?res
:res1
;
164 else return (res
>res1
)?res
:res1
;
167 /************************************************************************/
168 /*Routine receive double x and two double results of cos(x) and return */
169 /*result which is more accurate */
170 /*Computing cos(x) with multi precision routine c32 */
171 /************************************************************************/
174 __cos32(double x
, double res
, double res1
) {
179 __dbl_mp(0.5*(res1
-res
),&b
,p
);
182 { __sub(&pi
,&c
,&a
,p
);
187 { __sub(&hp
,&c
,&a
,p
);
190 else __c32(&c
,&b
,&a
,p
); /* b=cos(0.5*(res+res1)) */
191 __dbl_mp(x
,&c
,p
); /* c = x */
193 /* if a>0 return max(res,res1), otherwise return min(res,res1) */
194 if (a
.d
[0]>0) return (res
>res1
)?res
:res1
;
195 else return (res
<res1
)?res
:res1
;
198 /*******************************************************************/
199 /*Compute sin(x+dx) as Multi Precision number and return result as */
201 /*******************************************************************/
204 __mpsin(double x
, double dx
) {
212 if (x
>0.8) { __sub(&hp
,&c
,&a
,p
); __c32(&a
,&b
,&c
,p
); }
213 else __c32(&c
,&a
,&b
,p
); /* b = sin(x+dx) */
218 /*******************************************************************/
219 /* Compute cos()of double-length number (x+dx) as Multi Precision */
220 /* number and return result as double */
221 /*******************************************************************/
224 __mpcos(double x
, double dx
) {
233 { __sub(&hp
,&c
,&b
,p
);
236 else __c32(&c
,&a
,&b
,p
); /* a = cos(x+dx) */
241 /******************************************************************/
242 /* mpranred() performs range reduction of a double number x into */
243 /* multi precision number y, such that y=x-n*pi/2, abs(y)<pi/4, */
244 /* n=0,+-1,+-2,.... */
245 /* Return int which indicates in which quarter of circle x is */
246 /******************************************************************/
249 __mpranred(double x
, mp_no
*y
, int p
)
254 static const mp_no one
= {1,{1.0,1.0}};
257 if (ABS(x
) < 2.8e14
) {
258 t
= (x
*hpinv
.d
+ toint
.d
);
268 else { /* if x is very big more precision required */
275 for (i
=0;i
<p
;i
++) b
.d
[i
+1] = toverp
[i
+k
];
278 for (i
=1;i
<=p
-c
.e
;i
++) c
.d
[i
]=c
.d
[i
+c
.e
];
279 for (i
=p
+1-c
.e
;i
<=p
;i
++) c
.d
[i
]=0;
281 if (c
.d
[1] >= 8388608.0)
286 else __mul(&c
,&hp
,y
,p
);
288 if (x
< 0) { y
->d
[0] = - y
->d
[0]; n
= -n
; }
293 /*******************************************************************/
294 /* Multi-Precision sin() function subroutine, for p=32. It is */
295 /* based on the routines mpranred() and c32(). */
296 /*******************************************************************/
306 n
=__mpranred(x
,&u
,p
); /* n is 0, 1, 2 or 3 */
308 switch (n
) { /* in which quarter of unit circle y is*/
330 return 0; /* unreachable, to make the compiler happy */
333 /*****************************************************************/
334 /* Multi-Precision cos() function subroutine, for p=32. It is */
335 /* based on the routines mpranred() and c32(). */
336 /*****************************************************************/
348 n
=__mpranred(x
,&u
,p
); /* n is 0, 1, 2 or 3 */
350 switch (n
) { /* in what quarter of unit circle y is*/
373 return 0; /* unreachable, to make the compiler happy */
375 /******************************************************************/