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