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.  */
45 double
46 SECTION
__ieee754_asin(double x)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
libm_alias_finite(__ieee754_asin,__asin)178 libm_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.  */
188 double
189 SECTION
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
347 libm_alias_finite (__ieee754_acos, __acos)
348 #endif
349