1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2020 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 mpa.h mydefs.h usncs.h */ |
25 | /* doasin.c sincos32.c dosincos.c mpa.c */ |
26 | /* sincos.tbl asincos.tbl powtwo.tbl root.tbl */ |
27 | /* */ |
28 | /* Ultimate asin/acos routines. Given an IEEE double machine */ |
29 | /* number x, compute the correctly rounded value of */ |
30 | /* arcsin(x)or arccos(x) according to the function called. */ |
31 | /* Assumption: Machine arithmetic operations are performed in */ |
32 | /* round to nearest mode of IEEE 754 standard. */ |
33 | /* */ |
34 | /******************************************************************/ |
35 | #include "endian.h" |
36 | #include "mydefs.h" |
37 | #include "asincos.tbl" |
38 | #include "root.tbl" |
39 | #include "powtwo.tbl" |
40 | #include "MathLib.h" |
41 | #include "uasncs.h" |
42 | #include <float.h> |
43 | #include <math.h> |
44 | #include <math_private.h> |
45 | #include <math-underflow.h> |
46 | #include <libm-alias-finite.h> |
47 | |
48 | #ifndef SECTION |
49 | # define SECTION |
50 | #endif |
51 | |
52 | void __doasin(double x, double dx, double w[]); |
53 | void __dubsin(double x, double dx, double v[]); |
54 | void __dubcos(double x, double dx, double v[]); |
55 | void __docos(double x, double dx, double v[]); |
56 | double __sin32(double x, double res, double res1); |
57 | double __cos32(double x, double res, double res1); |
58 | |
59 | /***************************************************************************/ |
60 | /* An ultimate asin routine. Given an IEEE double machine number x */ |
61 | /* it computes the correctly rounded (to nearest) value of arcsin(x) */ |
62 | /***************************************************************************/ |
63 | double |
64 | SECTION |
65 | __ieee754_asin(double x){ |
66 | double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2]; |
67 | mynumber u,v; |
68 | int4 k,m,n; |
69 | |
70 | u.x = x; |
71 | m = u.i[HIGH_HALF]; |
72 | k = 0x7fffffff&m; /* no sign */ |
73 | |
74 | if (k < 0x3e500000) |
75 | { |
76 | math_check_force_underflow (x); |
77 | return x; /* for x->0 => sin(x)=x */ |
78 | } |
79 | /*----------------------2^-26 <= |x| < 2^ -3 -----------------*/ |
80 | else |
81 | if (k < 0x3fc00000) { |
82 | x2 = x*x; |
83 | t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x); |
84 | res = x+t; /* res=arcsin(x) according to Taylor series */ |
85 | cor = (x-res)+t; |
86 | if (res == res+1.025*cor) return res; |
87 | else { |
88 | x1 = x+big; |
89 | xx = x*x; |
90 | x1 -= big; |
91 | x2 = x - x1; |
92 | p = x1*x1*x1; |
93 | s1 = a1.x*p; |
94 | s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x + |
95 | ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p; |
96 | res1 = x+s1; |
97 | s2 = ((x-res1)+s1)+s2; |
98 | res = res1+s2; |
99 | cor = (res1-res)+s2; |
100 | if (res == res+1.00014*cor) return res; |
101 | else { |
102 | __doasin(x,0,w); |
103 | if (w[0]==(w[0]+1.00000001*w[1])) return w[0]; |
104 | else { |
105 | y=fabs(x); |
106 | res=fabs(w[0]); |
107 | res1=fabs(w[0]+1.1*w[1]); |
108 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
109 | } |
110 | } |
111 | } |
112 | } |
113 | /*---------------------0.125 <= |x| < 0.5 -----------------------------*/ |
114 | else if (k < 0x3fe00000) { |
115 | if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15); |
116 | else n = 11*((k&0x000fffff)>>14)+352; |
117 | if (m>0) xx = x - asncs.x[n]; |
118 | else xx = -x - asncs.x[n]; |
119 | t = asncs.x[n+1]*xx; |
120 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5] |
121 | +xx*asncs.x[n+6]))))+asncs.x[n+7]; |
122 | t+=p; |
123 | res =asncs.x[n+8] +t; |
124 | cor = (asncs.x[n+8]-res)+t; |
125 | if (res == res+1.05*cor) return (m>0)?res:-res; |
126 | else { |
127 | r=asncs.x[n+8]+xx*asncs.x[n+9]; |
128 | t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]); |
129 | res = r+t; |
130 | cor = (r-res)+t; |
131 | if (res == res+1.0005*cor) return (m>0)?res:-res; |
132 | else { |
133 | res1=res+1.1*cor; |
134 | z=0.5*(res1-res); |
135 | __dubsin(res,z,w); |
136 | z=(w[0]-fabs(x))+w[1]; |
137 | if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); |
138 | else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); |
139 | else { |
140 | y=fabs(x); |
141 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
142 | } |
143 | } |
144 | } |
145 | } /* else if (k < 0x3fe00000) */ |
146 | /*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/ |
147 | else |
148 | if (k < 0x3fe80000) { |
149 | n = 1056+((k&0x000fe000)>>11)*3; |
150 | if (m>0) xx = x - asncs.x[n]; |
151 | else xx = -x - asncs.x[n]; |
152 | t = asncs.x[n+1]*xx; |
153 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5] |
154 | +xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8]; |
155 | t+=p; |
156 | res =asncs.x[n+9] +t; |
157 | cor = (asncs.x[n+9]-res)+t; |
158 | if (res == res+1.01*cor) return (m>0)?res:-res; |
159 | else { |
160 | r=asncs.x[n+9]+xx*asncs.x[n+10]; |
161 | t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]); |
162 | res = r+t; |
163 | cor = (r-res)+t; |
164 | if (res == res+1.0005*cor) return (m>0)?res:-res; |
165 | else { |
166 | res1=res+1.1*cor; |
167 | z=0.5*(res1-res); |
168 | __dubsin(res,z,w); |
169 | z=(w[0]-fabs(x))+w[1]; |
170 | if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); |
171 | else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); |
172 | else { |
173 | y=fabs(x); |
174 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
175 | } |
176 | } |
177 | } |
178 | } /* else if (k < 0x3fe80000) */ |
179 | /*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/ |
180 | else |
181 | if (k < 0x3fed8000) { |
182 | n = 992+((k&0x000fe000)>>13)*13; |
183 | if (m>0) xx = x - asncs.x[n]; |
184 | else xx = -x - asncs.x[n]; |
185 | t = asncs.x[n+1]*xx; |
186 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5] |
187 | +xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9]; |
188 | t+=p; |
189 | res =asncs.x[n+10] +t; |
190 | cor = (asncs.x[n+10]-res)+t; |
191 | if (res == res+1.01*cor) return (m>0)?res:-res; |
192 | else { |
193 | r=asncs.x[n+10]+xx*asncs.x[n+11]; |
194 | t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]); |
195 | res = r+t; |
196 | cor = (r-res)+t; |
197 | if (res == res+1.0008*cor) return (m>0)?res:-res; |
198 | else { |
199 | res1=res+1.1*cor; |
200 | z=0.5*(res1-res); |
201 | y=hp0.x-res; |
202 | z=((hp0.x-y)-res)+(hp1.x-z); |
203 | __dubcos(y,z,w); |
204 | z=(w[0]-fabs(x))+w[1]; |
205 | if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); |
206 | else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); |
207 | else { |
208 | y=fabs(x); |
209 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
210 | } |
211 | } |
212 | } |
213 | } /* else if (k < 0x3fed8000) */ |
214 | /*-------------------0.921875 <= |x| < 0.953125 ------------------------*/ |
215 | else |
216 | if (k < 0x3fee8000) { |
217 | n = 884+((k&0x000fe000)>>13)*14; |
218 | if (m>0) xx = x - asncs.x[n]; |
219 | else xx = -x - asncs.x[n]; |
220 | t = asncs.x[n+1]*xx; |
221 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
222 | xx*(asncs.x[n+5]+xx*(asncs.x[n+6] |
223 | +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+ |
224 | xx*asncs.x[n+9])))))))+asncs.x[n+10]; |
225 | t+=p; |
226 | res =asncs.x[n+11] +t; |
227 | cor = (asncs.x[n+11]-res)+t; |
228 | if (res == res+1.01*cor) return (m>0)?res:-res; |
229 | else { |
230 | r=asncs.x[n+11]+xx*asncs.x[n+12]; |
231 | t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]); |
232 | res = r+t; |
233 | cor = (r-res)+t; |
234 | if (res == res+1.0007*cor) return (m>0)?res:-res; |
235 | else { |
236 | res1=res+1.1*cor; |
237 | z=0.5*(res1-res); |
238 | y=(hp0.x-res)-z; |
239 | z=y+hp1.x; |
240 | y=(y-z)+hp1.x; |
241 | __dubcos(z,y,w); |
242 | z=(w[0]-fabs(x))+w[1]; |
243 | if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); |
244 | else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); |
245 | else { |
246 | y=fabs(x); |
247 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
248 | } |
249 | } |
250 | } |
251 | } /* else if (k < 0x3fee8000) */ |
252 | |
253 | /*--------------------0.953125 <= |x| < 0.96875 ------------------------*/ |
254 | else |
255 | if (k < 0x3fef0000) { |
256 | n = 768+((k&0x000fe000)>>13)*15; |
257 | if (m>0) xx = x - asncs.x[n]; |
258 | else xx = -x - asncs.x[n]; |
259 | t = asncs.x[n+1]*xx; |
260 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
261 | xx*(asncs.x[n+5]+xx*(asncs.x[n+6] |
262 | +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+ |
263 | xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11]; |
264 | t+=p; |
265 | res =asncs.x[n+12] +t; |
266 | cor = (asncs.x[n+12]-res)+t; |
267 | if (res == res+1.01*cor) return (m>0)?res:-res; |
268 | else { |
269 | r=asncs.x[n+12]+xx*asncs.x[n+13]; |
270 | t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]); |
271 | res = r+t; |
272 | cor = (r-res)+t; |
273 | if (res == res+1.0007*cor) return (m>0)?res:-res; |
274 | else { |
275 | res1=res+1.1*cor; |
276 | z=0.5*(res1-res); |
277 | y=(hp0.x-res)-z; |
278 | z=y+hp1.x; |
279 | y=(y-z)+hp1.x; |
280 | __dubcos(z,y,w); |
281 | z=(w[0]-fabs(x))+w[1]; |
282 | if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1); |
283 | else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1); |
284 | else { |
285 | y=fabs(x); |
286 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
287 | } |
288 | } |
289 | } |
290 | } /* else if (k < 0x3fef0000) */ |
291 | /*--------------------0.96875 <= |x| < 1 --------------------------------*/ |
292 | else |
293 | if (k<0x3ff00000) { |
294 | z = 0.5*((m>0)?(1.0-x):(1.0+x)); |
295 | v.x=z; |
296 | k=v.i[HIGH_HALF]; |
297 | t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)]; |
298 | r=1.0-t*t*z; |
299 | t = t*(rt0+r*(rt1+r*(rt2+r*rt3))); |
300 | c=t*z; |
301 | t=c*(1.5-0.5*t*c); |
302 | y=(c+t24)-t24; |
303 | cc = (z-y*y)/(t+y); |
304 | p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z; |
305 | cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p; |
306 | res1 = hp0.x - 2.0*y; |
307 | res =res1 + cor; |
308 | if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res; |
309 | else { |
310 | c=y+cc; |
311 | cc=(y-c)+cc; |
312 | __doasin(c,cc,w); |
313 | res1=hp0.x-2.0*w[0]; |
314 | cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]); |
315 | res = res1+cor; |
316 | cor = (res1-res)+cor; |
317 | if (res==(res+1.0000001*cor)) return (m>0)?res:-res; |
318 | else { |
319 | y=fabs(x); |
320 | res1=res+1.1*cor; |
321 | return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1); |
322 | } |
323 | } |
324 | } /* else if (k < 0x3ff00000) */ |
325 | /*---------------------------- |x|>=1 -------------------------------*/ |
326 | else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x; |
327 | else |
328 | if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x; |
329 | else { |
330 | u.i[HIGH_HALF]=0x7ff00000; |
331 | v.i[HIGH_HALF]=0x7ff00000; |
332 | u.i[LOW_HALF]=0; |
333 | v.i[LOW_HALF]=0; |
334 | return u.x/v.x; /* NaN */ |
335 | } |
336 | } |
337 | #ifndef __ieee754_asin |
338 | libm_alias_finite (__ieee754_asin, __asin) |
339 | #endif |
340 | |
341 | /*******************************************************************/ |
342 | /* */ |
343 | /* End of arcsine, below is arccosine */ |
344 | /* */ |
345 | /*******************************************************************/ |
346 | |
347 | double |
348 | SECTION |
349 | __ieee754_acos(double x) |
350 | { |
351 | double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps; |
352 | mynumber u,v; |
353 | int4 k,m,n; |
354 | u.x = x; |
355 | m = u.i[HIGH_HALF]; |
356 | k = 0x7fffffff&m; |
357 | /*------------------- |x|<2.77556*10^-17 ----------------------*/ |
358 | if (k < 0x3c880000) return hp0.x; |
359 | |
360 | /*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/ |
361 | else |
362 | if (k < 0x3fc00000) { |
363 | x2 = x*x; |
364 | t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x); |
365 | r=hp0.x-x; |
366 | cor=(((hp0.x-r)-x)+hp1.x)-t; |
367 | res = r+cor; |
368 | cor = (r-res)+cor; |
369 | if (res == res+1.004*cor) return res; |
370 | else { |
371 | x1 = x+big; |
372 | xx = x*x; |
373 | x1 -= big; |
374 | x2 = x - x1; |
375 | p = x1*x1*x1; |
376 | s1 = a1.x*p; |
377 | s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x + |
378 | ((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p; |
379 | res1 = x+s1; |
380 | s2 = ((x-res1)+s1)+s2; |
381 | r=hp0.x-res1; |
382 | cor=(((hp0.x-r)-res1)+hp1.x)-s2; |
383 | res = r+cor; |
384 | cor = (r-res)+cor; |
385 | if (res == res+1.00004*cor) return res; |
386 | else { |
387 | __doasin(x,0,w); |
388 | r=hp0.x-w[0]; |
389 | cor=((hp0.x-r)-w[0])+(hp1.x-w[1]); |
390 | res=r+cor; |
391 | cor=(r-res)+cor; |
392 | if (res ==(res +1.00000001*cor)) return res; |
393 | else { |
394 | res1=res+1.1*cor; |
395 | return __cos32(x,res,res1); |
396 | } |
397 | } |
398 | } |
399 | } /* else if (k < 0x3fc00000) */ |
400 | /*---------------------- 0.125 <= |x| < 0.5 --------------------*/ |
401 | else |
402 | if (k < 0x3fe00000) { |
403 | if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15); |
404 | else n = 11*((k&0x000fffff)>>14)+352; |
405 | if (m>0) xx = x - asncs.x[n]; |
406 | else xx = -x - asncs.x[n]; |
407 | t = asncs.x[n+1]*xx; |
408 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
409 | xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7]; |
410 | t+=p; |
411 | y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]); |
412 | t = (m>0)?(hp1.x-t):(hp1.x+t); |
413 | res = y+t; |
414 | if (res == res+1.02*((y-res)+t)) return res; |
415 | else { |
416 | r=asncs.x[n+8]+xx*asncs.x[n+9]; |
417 | t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]); |
418 | if (m>0) |
419 | {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; } |
420 | else |
421 | {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); } |
422 | res = p+t; |
423 | cor = (p-res)+t; |
424 | if (res == (res+1.0002*cor)) return res; |
425 | else { |
426 | res1=res+1.1*cor; |
427 | z=0.5*(res1-res); |
428 | __docos(res,z,w); |
429 | z=(w[0]-x)+w[1]; |
430 | if (z>1.0e-27) return max(res,res1); |
431 | else if (z<-1.0e-27) return min(res,res1); |
432 | else return __cos32(x,res,res1); |
433 | } |
434 | } |
435 | } /* else if (k < 0x3fe00000) */ |
436 | |
437 | /*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/ |
438 | else |
439 | if (k < 0x3fe80000) { |
440 | n = 1056+((k&0x000fe000)>>11)*3; |
441 | if (m>0) {xx = x - asncs.x[n]; eps=1.04; } |
442 | else {xx = -x - asncs.x[n]; eps=1.02; } |
443 | t = asncs.x[n+1]*xx; |
444 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
445 | xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+ |
446 | xx*asncs.x[n+7])))))+asncs.x[n+8]; |
447 | t+=p; |
448 | y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]); |
449 | t = (m>0)?(hp1.x-t):(hp1.x+t); |
450 | res = y+t; |
451 | if (res == res+eps*((y-res)+t)) return res; |
452 | else { |
453 | r=asncs.x[n+9]+xx*asncs.x[n+10]; |
454 | t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]); |
455 | if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; } |
456 | else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; } |
457 | res = p+t; |
458 | cor = (p-res)+t; |
459 | if (res == (res+eps*cor)) return res; |
460 | else { |
461 | res1=res+1.1*cor; |
462 | z=0.5*(res1-res); |
463 | __docos(res,z,w); |
464 | z=(w[0]-x)+w[1]; |
465 | if (z>1.0e-27) return max(res,res1); |
466 | else if (z<-1.0e-27) return min(res,res1); |
467 | else return __cos32(x,res,res1); |
468 | } |
469 | } |
470 | } /* else if (k < 0x3fe80000) */ |
471 | |
472 | /*------------------------- 0.75 <= |x| < 0.921875 -------------*/ |
473 | else |
474 | if (k < 0x3fed8000) { |
475 | n = 992+((k&0x000fe000)>>13)*13; |
476 | if (m>0) {xx = x - asncs.x[n]; eps = 1.04; } |
477 | else {xx = -x - asncs.x[n]; eps = 1.01; } |
478 | t = asncs.x[n+1]*xx; |
479 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
480 | xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+ |
481 | xx*asncs.x[n+8]))))))+asncs.x[n+9]; |
482 | t+=p; |
483 | y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]); |
484 | t = (m>0)?(hp1.x-t):(hp1.x+t); |
485 | res = y+t; |
486 | if (res == res+eps*((y-res)+t)) return res; |
487 | else { |
488 | r=asncs.x[n+10]+xx*asncs.x[n+11]; |
489 | t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]); |
490 | if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; } |
491 | else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; } |
492 | res = p+t; |
493 | cor = (p-res)+t; |
494 | if (res == (res+eps*cor)) return res; |
495 | else { |
496 | res1=res+1.1*cor; |
497 | z=0.5*(res1-res); |
498 | __docos(res,z,w); |
499 | z=(w[0]-x)+w[1]; |
500 | if (z>1.0e-27) return max(res,res1); |
501 | else if (z<-1.0e-27) return min(res,res1); |
502 | else return __cos32(x,res,res1); |
503 | } |
504 | } |
505 | } /* else if (k < 0x3fed8000) */ |
506 | |
507 | /*-------------------0.921875 <= |x| < 0.953125 ------------------*/ |
508 | else |
509 | if (k < 0x3fee8000) { |
510 | n = 884+((k&0x000fe000)>>13)*14; |
511 | if (m>0) {xx = x - asncs.x[n]; eps=1.04; } |
512 | else {xx = -x - asncs.x[n]; eps =1.005; } |
513 | t = asncs.x[n+1]*xx; |
514 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
515 | xx*(asncs.x[n+5]+xx*(asncs.x[n+6] |
516 | +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+ |
517 | xx*asncs.x[n+9])))))))+asncs.x[n+10]; |
518 | t+=p; |
519 | y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]); |
520 | t = (m>0)?(hp1.x-t):(hp1.x+t); |
521 | res = y+t; |
522 | if (res == res+eps*((y-res)+t)) return res; |
523 | else { |
524 | r=asncs.x[n+11]+xx*asncs.x[n+12]; |
525 | t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]); |
526 | if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; } |
527 | else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; } |
528 | res = p+t; |
529 | cor = (p-res)+t; |
530 | if (res == (res+eps*cor)) return res; |
531 | else { |
532 | res1=res+1.1*cor; |
533 | z=0.5*(res1-res); |
534 | __docos(res,z,w); |
535 | z=(w[0]-x)+w[1]; |
536 | if (z>1.0e-27) return max(res,res1); |
537 | else if (z<-1.0e-27) return min(res,res1); |
538 | else return __cos32(x,res,res1); |
539 | } |
540 | } |
541 | } /* else if (k < 0x3fee8000) */ |
542 | |
543 | /*--------------------0.953125 <= |x| < 0.96875 ----------------*/ |
544 | else |
545 | if (k < 0x3fef0000) { |
546 | n = 768+((k&0x000fe000)>>13)*15; |
547 | if (m>0) {xx = x - asncs.x[n]; eps=1.04; } |
548 | else {xx = -x - asncs.x[n]; eps=1.005;} |
549 | t = asncs.x[n+1]*xx; |
550 | p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+ |
551 | xx*(asncs.x[n+5]+xx*(asncs.x[n+6] |
552 | +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+ |
553 | xx*asncs.x[n+10]))))))))+asncs.x[n+11]; |
554 | t+=p; |
555 | y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]); |
556 | t = (m>0)?(hp1.x-t):(hp1.x+t); |
557 | res = y+t; |
558 | if (res == res+eps*((y-res)+t)) return res; |
559 | else { |
560 | r=asncs.x[n+12]+xx*asncs.x[n+13]; |
561 | t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]); |
562 | if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; } |
563 | else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; } |
564 | res = p+t; |
565 | cor = (p-res)+t; |
566 | if (res == (res+eps*cor)) return res; |
567 | else { |
568 | res1=res+1.1*cor; |
569 | z=0.5*(res1-res); |
570 | __docos(res,z,w); |
571 | z=(w[0]-x)+w[1]; |
572 | if (z>1.0e-27) return max(res,res1); |
573 | else if (z<-1.0e-27) return min(res,res1); |
574 | else return __cos32(x,res,res1); |
575 | } |
576 | } |
577 | } /* else if (k < 0x3fef0000) */ |
578 | /*-----------------0.96875 <= |x| < 1 ---------------------------*/ |
579 | |
580 | else |
581 | if (k<0x3ff00000) { |
582 | z = 0.5*((m>0)?(1.0-x):(1.0+x)); |
583 | v.x=z; |
584 | k=v.i[HIGH_HALF]; |
585 | t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)]; |
586 | r=1.0-t*t*z; |
587 | t = t*(rt0+r*(rt1+r*(rt2+r*rt3))); |
588 | c=t*z; |
589 | t=c*(1.5-0.5*t*c); |
590 | y = (t27*c+c)-t27*c; |
591 | cc = (z-y*y)/(t+y); |
592 | p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z; |
593 | if (m<0) { |
594 | cor = (hp1.x - cc)-(y+cc)*p; |
595 | res1 = hp0.x - y; |
596 | res =res1 + cor; |
597 | if (res == res+1.002*((res1-res)+cor)) return (res+res); |
598 | else { |
599 | c=y+cc; |
600 | cc=(y-c)+cc; |
601 | __doasin(c,cc,w); |
602 | res1=hp0.x-w[0]; |
603 | cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]); |
604 | res = res1+cor; |
605 | cor = (res1-res)+cor; |
606 | if (res==(res+1.000001*cor)) return (res+res); |
607 | else { |
608 | res=res+res; |
609 | res1=res+1.2*cor; |
610 | return __cos32(x,res,res1); |
611 | } |
612 | } |
613 | } |
614 | else { |
615 | cor = cc+p*(y+cc); |
616 | res = y + cor; |
617 | if (res == res+1.03*((y-res)+cor)) return (res+res); |
618 | else { |
619 | c=y+cc; |
620 | cc=(y-c)+cc; |
621 | __doasin(c,cc,w); |
622 | res = w[0]; |
623 | cor=w[1]; |
624 | if (res==(res+1.000001*cor)) return (res+res); |
625 | else { |
626 | res=res+res; |
627 | res1=res+1.2*cor; |
628 | return __cos32(x,res,res1); |
629 | } |
630 | } |
631 | } |
632 | } /* else if (k < 0x3ff00000) */ |
633 | |
634 | /*---------------------------- |x|>=1 -----------------------*/ |
635 | else |
636 | if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x; |
637 | else |
638 | if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x; |
639 | else { |
640 | u.i[HIGH_HALF]=0x7ff00000; |
641 | v.i[HIGH_HALF]=0x7ff00000; |
642 | u.i[LOW_HALF]=0; |
643 | v.i[LOW_HALF]=0; |
644 | return u.x/v.x; |
645 | } |
646 | } |
647 | #ifndef __ieee754_acos |
648 | libm_alias_finite (__ieee754_acos, __acos) |
649 | #endif |
650 | |