2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2023 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 <https://www.gnu.org/licenses/>.
19 /******************************************************************/
20 /* MODULE_NAME:uasncs.c */
22 /* FUNCTIONS: uasin */
24 /* FILES NEEDED: dla.h endian.h mydefs.h usncs.h */
25 /* sincos.tbl asincos.tbl powtwo.tbl root.tbl */
27 /******************************************************************/
30 #include "asincos.tbl"
36 #include <math_private.h>
37 #include <math-underflow.h>
38 #include <libm-alias-finite.h>
44 /* asin with max ULP of ~0.516 based on random sampling. */
47 __ieee754_asin(double x
){
48 double x2
,xx
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
;
54 k
= 0x7fffffff&m
; /* no sign */
58 math_check_force_underflow (x
);
59 return x
; /* for x->0 => sin(x)=x */
61 /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
65 t
= (((((f6
*x2
+ f5
)*x2
+ f4
)*x2
+ f3
)*x2
+ f2
)*x2
+ f1
)*(x2
*x
);
66 res
= x
+t
; /* res=arcsin(x) according to Taylor series */
67 /* Max ULP is 0.513. */
70 /*---------------------0.125 <= |x| < 0.5 -----------------------------*/
71 else if (k
< 0x3fe00000) {
72 if (k
<0x3fd00000) n
= 11*((k
&0x000fffff)>>15);
73 else n
= 11*((k
&0x000fffff)>>14)+352;
74 if (m
>0) xx
= x
- asncs
.x
[n
];
75 else xx
= -x
- asncs
.x
[n
];
77 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
78 +xx
*asncs
.x
[n
+6]))))+asncs
.x
[n
+7];
81 /* Max ULP is 0.524. */
82 return (m
>0)?res
:-res
;
83 } /* else if (k < 0x3fe00000) */
84 /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
87 n
= 1056+((k
&0x000fe000)>>11)*3;
88 if (m
>0) xx
= x
- asncs
.x
[n
];
89 else xx
= -x
- asncs
.x
[n
];
91 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
92 +xx
*(asncs
.x
[n
+6]+xx
*asncs
.x
[n
+7])))))+asncs
.x
[n
+8];
95 /* Max ULP is 0.505. */
96 return (m
>0)?res
:-res
;
97 } /* else if (k < 0x3fe80000) */
98 /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
100 if (k
< 0x3fed8000) {
101 n
= 992+((k
&0x000fe000)>>13)*13;
102 if (m
>0) xx
= x
- asncs
.x
[n
];
103 else xx
= -x
- asncs
.x
[n
];
105 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+xx
*(asncs
.x
[n
+5]
106 +xx
*(asncs
.x
[n
+6]+xx
*(asncs
.x
[n
+7]+xx
*asncs
.x
[n
+8]))))))+asncs
.x
[n
+9];
108 res
=asncs
.x
[n
+10] +t
;
109 /* Max ULP is 0.505. */
110 return (m
>0)?res
:-res
;
111 } /* else if (k < 0x3fed8000) */
112 /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
114 if (k
< 0x3fee8000) {
115 n
= 884+((k
&0x000fe000)>>13)*14;
116 if (m
>0) xx
= x
- asncs
.x
[n
];
117 else xx
= -x
- asncs
.x
[n
];
119 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
120 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
121 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
122 xx
*asncs
.x
[n
+9])))))))+asncs
.x
[n
+10];
124 res
=asncs
.x
[n
+11] +t
;
125 /* Max ULP is 0.505. */
126 return (m
>0)?res
:-res
;
127 } /* else if (k < 0x3fee8000) */
129 /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
131 if (k
< 0x3fef0000) {
132 n
= 768+((k
&0x000fe000)>>13)*15;
133 if (m
>0) xx
= x
- asncs
.x
[n
];
134 else xx
= -x
- asncs
.x
[n
];
136 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
137 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
138 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
139 xx
*(asncs
.x
[n
+9]+xx
*asncs
.x
[n
+10]))))))))+asncs
.x
[n
+11];
141 res
=asncs
.x
[n
+12] +t
;
142 /* Max ULP is 0.505. */
143 return (m
>0)?res
:-res
;
144 } /* else if (k < 0x3fef0000) */
145 /*--------------------0.96875 <= |x| < 1 --------------------------------*/
148 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
151 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
153 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
158 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
159 cor
= (hp1
.x
- 2.0*cc
)-2.0*(y
+cc
)*p
;
160 res1
= hp0
.x
- 2.0*y
;
162 /* Max ULP is 0.5015. */
163 return (m
>0)?res
:-res
;
164 } /* else if (k < 0x3ff00000) */
165 /*---------------------------- |x|>=1 -------------------------------*/
166 else if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?hp0
.x
:-hp0
.x
;
168 return (x
- x
) / (x
- x
);
170 #ifndef __ieee754_asin
171 libm_alias_finite (__ieee754_asin
, __asin
)
174 /*******************************************************************/
176 /* End of arcsine, below is arccosine */
178 /*******************************************************************/
180 /* acos with max ULP of ~0.523 based on random sampling. */
183 __ieee754_acos(double x
)
185 double x2
,xx
,res1
,p
,t
,res
,r
,cor
,cc
,y
,c
,z
;
191 /*------------------- |x|<2.77556*10^-17 ----------------------*/
192 if (k
< 0x3c880000) return hp0
.x
;
194 /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
196 if (k
< 0x3fc00000) {
198 t
= (((((f6
*x2
+ f5
)*x2
+ f4
)*x2
+ f3
)*x2
+ f2
)*x2
+ f1
)*(x2
*x
);
200 cor
=(((hp0
.x
-r
)-x
)+hp1
.x
)-t
;
202 /* Max ULP is 0.502. */
204 } /* else if (k < 0x3fc00000) */
205 /*---------------------- 0.125 <= |x| < 0.5 --------------------*/
207 if (k
< 0x3fe00000) {
208 if (k
<0x3fd00000) n
= 11*((k
&0x000fffff)>>15);
209 else n
= 11*((k
&0x000fffff)>>14)+352;
210 if (m
>0) xx
= x
- asncs
.x
[n
];
211 else xx
= -x
- asncs
.x
[n
];
213 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
214 xx
*(asncs
.x
[n
+5]+xx
*asncs
.x
[n
+6]))))+asncs
.x
[n
+7];
216 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+8]):(hp0
.x
+asncs
.x
[n
+8]);
217 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
219 /* Max ULP is 0.51. */
221 } /* else if (k < 0x3fe00000) */
223 /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
225 if (k
< 0x3fe80000) {
226 n
= 1056+((k
&0x000fe000)>>11)*3;
227 if (m
>0) {xx
= x
- asncs
.x
[n
]; }
228 else {xx
= -x
- asncs
.x
[n
]; }
230 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
231 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]+
232 xx
*asncs
.x
[n
+7])))))+asncs
.x
[n
+8];
234 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+9]):(hp0
.x
+asncs
.x
[n
+9]);
235 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
237 /* Max ULP is 0.523 based on random sampling. */
239 } /* else if (k < 0x3fe80000) */
241 /*------------------------- 0.75 <= |x| < 0.921875 -------------*/
243 if (k
< 0x3fed8000) {
244 n
= 992+((k
&0x000fe000)>>13)*13;
245 if (m
>0) {xx
= x
- asncs
.x
[n
]; }
246 else {xx
= -x
- asncs
.x
[n
]; }
248 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
249 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]+xx
*(asncs
.x
[n
+7]+
250 xx
*asncs
.x
[n
+8]))))))+asncs
.x
[n
+9];
252 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+10]):(hp0
.x
+asncs
.x
[n
+10]);
253 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
255 /* Max ULP is 0.523 based on random sampling. */
257 } /* else if (k < 0x3fed8000) */
259 /*-------------------0.921875 <= |x| < 0.953125 ------------------*/
261 if (k
< 0x3fee8000) {
262 n
= 884+((k
&0x000fe000)>>13)*14;
263 if (m
>0) {xx
= x
- asncs
.x
[n
]; }
264 else {xx
= -x
- asncs
.x
[n
]; }
266 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
267 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
268 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+
269 xx
*asncs
.x
[n
+9])))))))+asncs
.x
[n
+10];
271 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+11]):(hp0
.x
+asncs
.x
[n
+11]);
272 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
274 /* Max ULP is 0.523 based on random sampling. */
276 } /* else if (k < 0x3fee8000) */
278 /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
280 if (k
< 0x3fef0000) {
281 n
= 768+((k
&0x000fe000)>>13)*15;
282 if (m
>0) {xx
= x
- asncs
.x
[n
]; }
283 else {xx
= -x
- asncs
.x
[n
]; }
285 p
=xx
*xx
*(asncs
.x
[n
+2]+xx
*(asncs
.x
[n
+3]+xx
*(asncs
.x
[n
+4]+
286 xx
*(asncs
.x
[n
+5]+xx
*(asncs
.x
[n
+6]
287 +xx
*(asncs
.x
[n
+7]+xx
*(asncs
.x
[n
+8]+xx
*(asncs
.x
[n
+9]+
288 xx
*asncs
.x
[n
+10]))))))))+asncs
.x
[n
+11];
290 y
= (m
>0)?(hp0
.x
-asncs
.x
[n
+12]):(hp0
.x
+asncs
.x
[n
+12]);
291 t
= (m
>0)?(hp1
.x
-t
):(hp1
.x
+t
);
293 /* Max ULP is 0.523 based on random sampling. */
295 } /* else if (k < 0x3fef0000) */
296 /*-----------------0.96875 <= |x| < 1 ---------------------------*/
300 z
= 0.5*((m
>0)?(1.0-x
):(1.0+x
));
303 t
=inroot
[(k
&0x001fffff)>>14]*powtwo
[511-(k
>>21)];
305 t
= t
*(rt0
+r
*(rt1
+r
*(rt2
+r
*rt3
)));
310 p
=(((((f6
*z
+f5
)*z
+f4
)*z
+f3
)*z
+f2
)*z
+f1
)*z
;
312 cor
= (hp1
.x
- cc
)-(y
+cc
)*p
;
315 /* Max ULP is 0.501. */
321 /* Max ULP is 0.515. */
324 } /* else if (k < 0x3ff00000) */
326 /*---------------------------- |x|>=1 -----------------------*/
328 if (k
==0x3ff00000 && u
.i
[LOW_HALF
]==0) return (m
>0)?0:2.0*hp0
.x
;
330 return (x
- x
) / (x
- x
);
332 #ifndef __ieee754_acos
333 libm_alias_finite (__ieee754_acos
, __acos
)