1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2021 Free Software Foundation, Inc.
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
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.
15 *
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/>.
18 */
19/******************************************************************/
20/* MODULE_NAME:uasncs.c */
21/* */
22/* FUNCTIONS: uasin */
23/* uacos */
24/* FILES NEEDED: dla.h endian.h mydefs.h usncs.h */
25/* sincos.tbl asincos.tbl powtwo.tbl root.tbl */
26/* */
27/******************************************************************/
28#include "endian.h"
29#include "mydefs.h"
30#include "asincos.tbl"
31#include "root.tbl"
32#include "powtwo.tbl"
33#include "uasncs.h"
34#include <float.h>
35#include <math.h>
36#include <math_private.h>
37#include <math-underflow.h>
38#include <libm-alias-finite.h>
39
40#ifndef SECTION
41# define SECTION
42#endif
43
44/* asin with max ULP of ~0.516 based on random sampling. */
45double
46SECTION
47__ieee754_asin(double x){
48 double x2,xx,res1,p,t,res,r,cor,cc,y,c,z;
49 mynumber u,v;
50 int4 k,m,n;
51
52 u.x = x;
53 m = u.i[HIGH_HALF];
54 k = 0x7fffffff&m; /* no sign */
55
56 if (k < 0x3e500000)
57 {
58 math_check_force_underflow (x);
59 return x; /* for x->0 => sin(x)=x */
60 }
61 /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
62 else
63 if (k < 0x3fc00000) {
64 x2 = x*x;
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. */
68 return res;
69 }
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];
76 t = asncs.x[n+1]*xx;
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];
79 t+=p;
80 res =asncs.x[n+8] +t;
81 /* Max ULP is 0.524. */
82 return (m>0)?res:-res;
83 } /* else if (k < 0x3fe00000) */
84 /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
85 else
86 if (k < 0x3fe80000) {
87 n = 1056+((k&0x000fe000)>>11)*3;
88 if (m>0) xx = x - asncs.x[n];
89 else xx = -x - asncs.x[n];
90 t = asncs.x[n+1]*xx;
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];
93 t+=p;
94 res =asncs.x[n+9] +t;
95 /* Max ULP is 0.505. */
96 return (m>0)?res:-res;
97 } /* else if (k < 0x3fe80000) */
98 /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
99 else
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];
104 t = asncs.x[n+1]*xx;
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];
107 t+=p;
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 ------------------------*/
113 else
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];
118 t = asncs.x[n+1]*xx;
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];
123 t+=p;
124 res =asncs.x[n+11] +t;
125 /* Max ULP is 0.505. */
126 return (m>0)?res:-res;
127 } /* else if (k < 0x3fee8000) */
128
129 /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
130 else
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];
135 t = asncs.x[n+1]*xx;
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];
140 t+=p;
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 --------------------------------*/
146 else
147 if (k<0x3ff00000) {
148 z = 0.5*((m>0)?(1.0-x):(1.0+x));
149 v.x=z;
150 k=v.i[HIGH_HALF];
151 t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
152 r=1.0-t*t*z;
153 t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
154 c=t*z;
155 t=c*(1.5-0.5*t*c);
156 y=(c+t24)-t24;
157 cc = (z-y*y)/(t+y);
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;
161 res =res1 + cor;
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;
167 else
168 if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x;
169 else {
170 u.i[HIGH_HALF]=0x7ff00000;
171 v.i[HIGH_HALF]=0x7ff00000;
172 u.i[LOW_HALF]=0;
173 v.i[LOW_HALF]=0;
174 return u.x/v.x; /* NaN */
175 }
176}
177#ifndef __ieee754_asin
178libm_alias_finite (__ieee754_asin, __asin)
179#endif
180
181/*******************************************************************/
182/* */
183/* End of arcsine, below is arccosine */
184/* */
185/*******************************************************************/
186
187/* acos with max ULP of ~0.523 based on random sampling. */
188double
189SECTION
190__ieee754_acos(double x)
191{
192 double x2,xx,res1,p,t,res,r,cor,cc,y,c,z;
193 mynumber u,v;
194 int4 k,m,n;
195 u.x = x;
196 m = u.i[HIGH_HALF];
197 k = 0x7fffffff&m;
198 /*------------------- |x|<2.77556*10^-17 ----------------------*/
199 if (k < 0x3c880000) return hp0.x;
200
201 /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
202 else
203 if (k < 0x3fc00000) {
204 x2 = x*x;
205 t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
206 r=hp0.x-x;
207 cor=(((hp0.x-r)-x)+hp1.x)-t;
208 res = r+cor;
209 /* Max ULP is 0.502. */
210 return res;
211 } /* else if (k < 0x3fc00000) */
212 /*---------------------- 0.125 <= |x| < 0.5 --------------------*/
213 else
214 if (k < 0x3fe00000) {
215 if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
216 else n = 11*((k&0x000fffff)>>14)+352;
217 if (m>0) xx = x - asncs.x[n];
218 else xx = -x - asncs.x[n];
219 t = asncs.x[n+1]*xx;
220 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
221 xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
222 t+=p;
223 y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
224 t = (m>0)?(hp1.x-t):(hp1.x+t);
225 res = y+t;
226 /* Max ULP is 0.51. */
227 return res;
228 } /* else if (k < 0x3fe00000) */
229
230 /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
231 else
232 if (k < 0x3fe80000) {
233 n = 1056+((k&0x000fe000)>>11)*3;
234 if (m>0) {xx = x - asncs.x[n]; }
235 else {xx = -x - asncs.x[n]; }
236 t = asncs.x[n+1]*xx;
237 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
238 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
239 xx*asncs.x[n+7])))))+asncs.x[n+8];
240 t+=p;
241 y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
242 t = (m>0)?(hp1.x-t):(hp1.x+t);
243 res = y+t;
244 /* Max ULP is 0.523 based on random sampling. */
245 return res;
246 } /* else if (k < 0x3fe80000) */
247
248/*------------------------- 0.75 <= |x| < 0.921875 -------------*/
249 else
250 if (k < 0x3fed8000) {
251 n = 992+((k&0x000fe000)>>13)*13;
252 if (m>0) {xx = x - asncs.x[n]; }
253 else {xx = -x - asncs.x[n]; }
254 t = asncs.x[n+1]*xx;
255 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
256 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
257 xx*asncs.x[n+8]))))))+asncs.x[n+9];
258 t+=p;
259 y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
260 t = (m>0)?(hp1.x-t):(hp1.x+t);
261 res = y+t;
262 /* Max ULP is 0.523 based on random sampling. */
263 return res;
264 } /* else if (k < 0x3fed8000) */
265
266/*-------------------0.921875 <= |x| < 0.953125 ------------------*/
267 else
268 if (k < 0x3fee8000) {
269 n = 884+((k&0x000fe000)>>13)*14;
270 if (m>0) {xx = x - asncs.x[n]; }
271 else {xx = -x - asncs.x[n]; }
272 t = asncs.x[n+1]*xx;
273 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
274 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
275 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
276 xx*asncs.x[n+9])))))))+asncs.x[n+10];
277 t+=p;
278 y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
279 t = (m>0)?(hp1.x-t):(hp1.x+t);
280 res = y+t;
281 /* Max ULP is 0.523 based on random sampling. */
282 return res;
283 } /* else if (k < 0x3fee8000) */
284
285 /*--------------------0.953125 <= |x| < 0.96875 ----------------*/
286 else
287 if (k < 0x3fef0000) {
288 n = 768+((k&0x000fe000)>>13)*15;
289 if (m>0) {xx = x - asncs.x[n]; }
290 else {xx = -x - asncs.x[n]; }
291 t = asncs.x[n+1]*xx;
292 p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
293 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
294 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
295 xx*asncs.x[n+10]))))))))+asncs.x[n+11];
296 t+=p;
297 y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
298 t = (m>0)?(hp1.x-t):(hp1.x+t);
299 res = y+t;
300 /* Max ULP is 0.523 based on random sampling. */
301 return res;
302 } /* else if (k < 0x3fef0000) */
303 /*-----------------0.96875 <= |x| < 1 ---------------------------*/
304
305 else
306 if (k<0x3ff00000) {
307 z = 0.5*((m>0)?(1.0-x):(1.0+x));
308 v.x=z;
309 k=v.i[HIGH_HALF];
310 t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
311 r=1.0-t*t*z;
312 t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
313 c=t*z;
314 t=c*(1.5-0.5*t*c);
315 y = (t27*c+c)-t27*c;
316 cc = (z-y*y)/(t+y);
317 p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
318 if (m<0) {
319 cor = (hp1.x - cc)-(y+cc)*p;
320 res1 = hp0.x - y;
321 res =res1 + cor;
322 /* Max ULP is 0.501. */
323 return (res+res);
324 }
325 else {
326 cor = cc+p*(y+cc);
327 res = y + cor;
328 /* Max ULP is 0.515. */
329 return (res+res);
330 }
331 } /* else if (k < 0x3ff00000) */
332
333 /*---------------------------- |x|>=1 -----------------------*/
334 else
335 if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
336 else
337 if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x;
338 else {
339 u.i[HIGH_HALF]=0x7ff00000;
340 v.i[HIGH_HALF]=0x7ff00000;
341 u.i[LOW_HALF]=0;
342 v.i[LOW_HALF]=0;
343 return u.x/v.x;
344 }
345}
346#ifndef __ieee754_acos
347libm_alias_finite (__ieee754_acos, __acos)
348#endif
349