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