mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-30 08:40:07 +00:00
ec0ce0d3be
Similar to various other bugs in this area, some asin implementations do not raise the underflow exception for subnormal arguments, when the result is tiny and inexact. This patch forces the exception in a similar way to previous fixes. Tested for x86_64, x86, powerpc and mips64. [BZ #16351] * sysdeps/i386/fpu/e_asin.S (dbl_min): New object. (MO): New macro. (__ieee754_asin): Force underflow exception for results with small absolute value. * sysdeps/i386/fpu/e_asinf.S (flt_min): New object. (MO): New macro. (__ieee754_asinf): Force underflow exception for results with small absolute value. * sysdeps/ieee754/dbl-64/e_asin.c: Include <float.h> and <math.h>. (__ieee754_asin): Force underflow exception for results with small absolute value. * sysdeps/ieee754/flt-32/e_asinf.c: Include <float.h>. (__ieee754_asinf): Force underflow exception for results with small absolute value. * sysdeps/ieee754/ldbl-128/e_asinl.c: Include <float.h>. (__ieee754_asinl): Force underflow exception for results with small absolute value. * sysdeps/ieee754/ldbl-128ibm/e_asinl.c: Include <float.h>. (__ieee754_asinl): Force underflow exception for results with small absolute value. * sysdeps/ieee754/ldbl-96/e_asinl.c: Include <float.h>. (__ieee754_asinl): Force underflow exception for results with small absolute value. * sysdeps/x86_64/fpu/multiarch/e_asin.c [HAVE_FMA4_SUPPORT]: Include <math.h>. * math/auto-libm-test-in: Do not mark underflow exceptions as possibly missing for bug 16351. * math/auto-libm-test-out: Regenerated.
652 lines
20 KiB
C
652 lines
20 KiB
C
/*
|
|
* IBM Accurate Mathematical Library
|
|
* written by International Business Machines Corp.
|
|
* Copyright (C) 2001-2015 Free Software Foundation, Inc.
|
|
*
|
|
* This program is free software; you can redistribute it and/or modify
|
|
* it under the terms of the GNU Lesser General Public License as published by
|
|
* the Free Software Foundation; either version 2.1 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Lesser General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public License
|
|
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
/******************************************************************/
|
|
/* MODULE_NAME:uasncs.c */
|
|
/* */
|
|
/* FUNCTIONS: uasin */
|
|
/* uacos */
|
|
/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
|
|
/* doasin.c sincos32.c dosincos.c mpa.c */
|
|
/* sincos.tbl asincos.tbl powtwo.tbl root.tbl */
|
|
/* */
|
|
/* Ultimate asin/acos routines. Given an IEEE double machine */
|
|
/* number x, compute the correctly rounded value of */
|
|
/* arcsin(x)or arccos(x) according to the function called. */
|
|
/* Assumption: Machine arithmetic operations are performed in */
|
|
/* round to nearest mode of IEEE 754 standard. */
|
|
/* */
|
|
/******************************************************************/
|
|
#include "endian.h"
|
|
#include "mydefs.h"
|
|
#include "asincos.tbl"
|
|
#include "root.tbl"
|
|
#include "powtwo.tbl"
|
|
#include "MathLib.h"
|
|
#include "uasncs.h"
|
|
#include <float.h>
|
|
#include <math.h>
|
|
#include <math_private.h>
|
|
|
|
#ifndef SECTION
|
|
# define SECTION
|
|
#endif
|
|
|
|
void __doasin(double x, double dx, double w[]);
|
|
void __dubsin(double x, double dx, double v[]);
|
|
void __dubcos(double x, double dx, double v[]);
|
|
void __docos(double x, double dx, double v[]);
|
|
double __sin32(double x, double res, double res1);
|
|
double __cos32(double x, double res, double res1);
|
|
|
|
/***************************************************************************/
|
|
/* An ultimate asin routine. Given an IEEE double machine number x */
|
|
/* it computes the correctly rounded (to nearest) value of arcsin(x) */
|
|
/***************************************************************************/
|
|
double
|
|
SECTION
|
|
__ieee754_asin(double x){
|
|
double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
|
|
mynumber u,v;
|
|
int4 k,m,n;
|
|
|
|
u.x = x;
|
|
m = u.i[HIGH_HALF];
|
|
k = 0x7fffffff&m; /* no sign */
|
|
|
|
if (k < 0x3e500000)
|
|
{
|
|
if (fabs (x) < DBL_MIN)
|
|
{
|
|
double force_underflow = x * x;
|
|
math_force_eval (force_underflow);
|
|
}
|
|
return x; /* for x->0 => sin(x)=x */
|
|
}
|
|
/*----------------------2^-26 <= |x| < 2^ -3 -----------------*/
|
|
else
|
|
if (k < 0x3fc00000) {
|
|
x2 = x*x;
|
|
t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
|
|
res = x+t; /* res=arcsin(x) according to Taylor series */
|
|
cor = (x-res)+t;
|
|
if (res == res+1.025*cor) return res;
|
|
else {
|
|
x1 = x+big;
|
|
xx = x*x;
|
|
x1 -= big;
|
|
x2 = x - x1;
|
|
p = x1*x1*x1;
|
|
s1 = a1.x*p;
|
|
s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
|
|
((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
|
|
res1 = x+s1;
|
|
s2 = ((x-res1)+s1)+s2;
|
|
res = res1+s2;
|
|
cor = (res1-res)+s2;
|
|
if (res == res+1.00014*cor) return res;
|
|
else {
|
|
__doasin(x,0,w);
|
|
if (w[0]==(w[0]+1.00000001*w[1])) return w[0];
|
|
else {
|
|
y=ABS(x);
|
|
res=ABS(w[0]);
|
|
res1=ABS(w[0]+1.1*w[1]);
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*---------------------0.125 <= |x| < 0.5 -----------------------------*/
|
|
else if (k < 0x3fe00000) {
|
|
if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
|
|
else n = 11*((k&0x000fffff)>>14)+352;
|
|
if (m>0) xx = x - asncs.x[n];
|
|
else xx = -x - asncs.x[n];
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
|
|
+xx*asncs.x[n+6]))))+asncs.x[n+7];
|
|
t+=p;
|
|
res =asncs.x[n+8] +t;
|
|
cor = (asncs.x[n+8]-res)+t;
|
|
if (res == res+1.05*cor) return (m>0)?res:-res;
|
|
else {
|
|
r=asncs.x[n+8]+xx*asncs.x[n+9];
|
|
t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
|
|
res = r+t;
|
|
cor = (r-res)+t;
|
|
if (res == res+1.0005*cor) return (m>0)?res:-res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__dubsin(res,z,w);
|
|
z=(w[0]-ABS(x))+w[1];
|
|
if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
|
|
else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
|
|
else {
|
|
y=ABS(x);
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fe00000) */
|
|
/*-------------------- 0.5 <= |x| < 0.75 -----------------------------*/
|
|
else
|
|
if (k < 0x3fe80000) {
|
|
n = 1056+((k&0x000fe000)>>11)*3;
|
|
if (m>0) xx = x - asncs.x[n];
|
|
else xx = -x - asncs.x[n];
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
|
|
+xx*(asncs.x[n+6]+xx*asncs.x[n+7])))))+asncs.x[n+8];
|
|
t+=p;
|
|
res =asncs.x[n+9] +t;
|
|
cor = (asncs.x[n+9]-res)+t;
|
|
if (res == res+1.01*cor) return (m>0)?res:-res;
|
|
else {
|
|
r=asncs.x[n+9]+xx*asncs.x[n+10];
|
|
t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
|
|
res = r+t;
|
|
cor = (r-res)+t;
|
|
if (res == res+1.0005*cor) return (m>0)?res:-res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__dubsin(res,z,w);
|
|
z=(w[0]-ABS(x))+w[1];
|
|
if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
|
|
else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
|
|
else {
|
|
y=ABS(x);
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fe80000) */
|
|
/*--------------------- 0.75 <= |x|< 0.921875 ----------------------*/
|
|
else
|
|
if (k < 0x3fed8000) {
|
|
n = 992+((k&0x000fe000)>>13)*13;
|
|
if (m>0) xx = x - asncs.x[n];
|
|
else xx = -x - asncs.x[n];
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+xx*(asncs.x[n+5]
|
|
+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+xx*asncs.x[n+8]))))))+asncs.x[n+9];
|
|
t+=p;
|
|
res =asncs.x[n+10] +t;
|
|
cor = (asncs.x[n+10]-res)+t;
|
|
if (res == res+1.01*cor) return (m>0)?res:-res;
|
|
else {
|
|
r=asncs.x[n+10]+xx*asncs.x[n+11];
|
|
t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
|
|
res = r+t;
|
|
cor = (r-res)+t;
|
|
if (res == res+1.0008*cor) return (m>0)?res:-res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
y=hp0.x-res;
|
|
z=((hp0.x-y)-res)+(hp1.x-z);
|
|
__dubcos(y,z,w);
|
|
z=(w[0]-ABS(x))+w[1];
|
|
if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
|
|
else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
|
|
else {
|
|
y=ABS(x);
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fed8000) */
|
|
/*-------------------0.921875 <= |x| < 0.953125 ------------------------*/
|
|
else
|
|
if (k < 0x3fee8000) {
|
|
n = 884+((k&0x000fe000)>>13)*14;
|
|
if (m>0) xx = x - asncs.x[n];
|
|
else xx = -x - asncs.x[n];
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
|
|
+xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
|
|
xx*asncs.x[n+9])))))))+asncs.x[n+10];
|
|
t+=p;
|
|
res =asncs.x[n+11] +t;
|
|
cor = (asncs.x[n+11]-res)+t;
|
|
if (res == res+1.01*cor) return (m>0)?res:-res;
|
|
else {
|
|
r=asncs.x[n+11]+xx*asncs.x[n+12];
|
|
t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
|
|
res = r+t;
|
|
cor = (r-res)+t;
|
|
if (res == res+1.0007*cor) return (m>0)?res:-res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
y=(hp0.x-res)-z;
|
|
z=y+hp1.x;
|
|
y=(y-z)+hp1.x;
|
|
__dubcos(z,y,w);
|
|
z=(w[0]-ABS(x))+w[1];
|
|
if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
|
|
else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
|
|
else {
|
|
y=ABS(x);
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fee8000) */
|
|
|
|
/*--------------------0.953125 <= |x| < 0.96875 ------------------------*/
|
|
else
|
|
if (k < 0x3fef0000) {
|
|
n = 768+((k&0x000fe000)>>13)*15;
|
|
if (m>0) xx = x - asncs.x[n];
|
|
else xx = -x - asncs.x[n];
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
|
|
+xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
|
|
xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
|
|
t+=p;
|
|
res =asncs.x[n+12] +t;
|
|
cor = (asncs.x[n+12]-res)+t;
|
|
if (res == res+1.01*cor) return (m>0)?res:-res;
|
|
else {
|
|
r=asncs.x[n+12]+xx*asncs.x[n+13];
|
|
t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
|
|
res = r+t;
|
|
cor = (r-res)+t;
|
|
if (res == res+1.0007*cor) return (m>0)?res:-res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
y=(hp0.x-res)-z;
|
|
z=y+hp1.x;
|
|
y=(y-z)+hp1.x;
|
|
__dubcos(z,y,w);
|
|
z=(w[0]-ABS(x))+w[1];
|
|
if (z>1.0e-27) return (m>0)?min(res,res1):-min(res,res1);
|
|
else if (z<-1.0e-27) return (m>0)?max(res,res1):-max(res,res1);
|
|
else {
|
|
y=ABS(x);
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fef0000) */
|
|
/*--------------------0.96875 <= |x| < 1 --------------------------------*/
|
|
else
|
|
if (k<0x3ff00000) {
|
|
z = 0.5*((m>0)?(1.0-x):(1.0+x));
|
|
v.x=z;
|
|
k=v.i[HIGH_HALF];
|
|
t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
|
|
r=1.0-t*t*z;
|
|
t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
|
|
c=t*z;
|
|
t=c*(1.5-0.5*t*c);
|
|
y=(c+t24)-t24;
|
|
cc = (z-y*y)/(t+y);
|
|
p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
|
|
cor = (hp1.x - 2.0*cc)-2.0*(y+cc)*p;
|
|
res1 = hp0.x - 2.0*y;
|
|
res =res1 + cor;
|
|
if (res == res+1.003*((res1-res)+cor)) return (m>0)?res:-res;
|
|
else {
|
|
c=y+cc;
|
|
cc=(y-c)+cc;
|
|
__doasin(c,cc,w);
|
|
res1=hp0.x-2.0*w[0];
|
|
cor=((hp0.x-res1)-2.0*w[0])+(hp1.x-2.0*w[1]);
|
|
res = res1+cor;
|
|
cor = (res1-res)+cor;
|
|
if (res==(res+1.0000001*cor)) return (m>0)?res:-res;
|
|
else {
|
|
y=ABS(x);
|
|
res1=res+1.1*cor;
|
|
return (m>0)?__sin32(y,res,res1):-__sin32(y,res,res1);
|
|
}
|
|
}
|
|
} /* else if (k < 0x3ff00000) */
|
|
/*---------------------------- |x|>=1 -------------------------------*/
|
|
else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
|
|
else
|
|
if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
|
|
else {
|
|
u.i[HIGH_HALF]=0x7ff00000;
|
|
v.i[HIGH_HALF]=0x7ff00000;
|
|
u.i[LOW_HALF]=0;
|
|
v.i[LOW_HALF]=0;
|
|
return u.x/v.x; /* NaN */
|
|
}
|
|
}
|
|
#ifndef __ieee754_asin
|
|
strong_alias (__ieee754_asin, __asin_finite)
|
|
#endif
|
|
|
|
/*******************************************************************/
|
|
/* */
|
|
/* End of arcsine, below is arccosine */
|
|
/* */
|
|
/*******************************************************************/
|
|
|
|
double
|
|
SECTION
|
|
__ieee754_acos(double x)
|
|
{
|
|
double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
|
|
mynumber u,v;
|
|
int4 k,m,n;
|
|
u.x = x;
|
|
m = u.i[HIGH_HALF];
|
|
k = 0x7fffffff&m;
|
|
/*------------------- |x|<2.77556*10^-17 ----------------------*/
|
|
if (k < 0x3c880000) return hp0.x;
|
|
|
|
/*----------------- 2.77556*10^-17 <= |x| < 2^-3 --------------*/
|
|
else
|
|
if (k < 0x3fc00000) {
|
|
x2 = x*x;
|
|
t = (((((f6*x2 + f5)*x2 + f4)*x2 + f3)*x2 + f2)*x2 + f1)*(x2*x);
|
|
r=hp0.x-x;
|
|
cor=(((hp0.x-r)-x)+hp1.x)-t;
|
|
res = r+cor;
|
|
cor = (r-res)+cor;
|
|
if (res == res+1.004*cor) return res;
|
|
else {
|
|
x1 = x+big;
|
|
xx = x*x;
|
|
x1 -= big;
|
|
x2 = x - x1;
|
|
p = x1*x1*x1;
|
|
s1 = a1.x*p;
|
|
s2 = ((((((c7*xx + c6)*xx + c5)*xx + c4)*xx + c3)*xx + c2)*xx*xx*x +
|
|
((a1.x+a2.x)*x2*x2+ 0.5*x1*x)*x2) + a2.x*p;
|
|
res1 = x+s1;
|
|
s2 = ((x-res1)+s1)+s2;
|
|
r=hp0.x-res1;
|
|
cor=(((hp0.x-r)-res1)+hp1.x)-s2;
|
|
res = r+cor;
|
|
cor = (r-res)+cor;
|
|
if (res == res+1.00004*cor) return res;
|
|
else {
|
|
__doasin(x,0,w);
|
|
r=hp0.x-w[0];
|
|
cor=((hp0.x-r)-w[0])+(hp1.x-w[1]);
|
|
res=r+cor;
|
|
cor=(r-res)+cor;
|
|
if (res ==(res +1.00000001*cor)) return res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fc00000) */
|
|
/*---------------------- 0.125 <= |x| < 0.5 --------------------*/
|
|
else
|
|
if (k < 0x3fe00000) {
|
|
if (k<0x3fd00000) n = 11*((k&0x000fffff)>>15);
|
|
else n = 11*((k&0x000fffff)>>14)+352;
|
|
if (m>0) xx = x - asncs.x[n];
|
|
else xx = -x - asncs.x[n];
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
|
|
t+=p;
|
|
y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
|
|
t = (m>0)?(hp1.x-t):(hp1.x+t);
|
|
res = y+t;
|
|
if (res == res+1.02*((y-res)+t)) return res;
|
|
else {
|
|
r=asncs.x[n+8]+xx*asncs.x[n+9];
|
|
t=((asncs.x[n+8]-r)+xx*asncs.x[n+9])+(p+xx*asncs.x[n+10]);
|
|
if (m>0)
|
|
{p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; }
|
|
else
|
|
{p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); }
|
|
res = p+t;
|
|
cor = (p-res)+t;
|
|
if (res == (res+1.0002*cor)) return res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__docos(res,z,w);
|
|
z=(w[0]-x)+w[1];
|
|
if (z>1.0e-27) return max(res,res1);
|
|
else if (z<-1.0e-27) return min(res,res1);
|
|
else return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fe00000) */
|
|
|
|
/*--------------------------- 0.5 <= |x| < 0.75 ---------------------*/
|
|
else
|
|
if (k < 0x3fe80000) {
|
|
n = 1056+((k&0x000fe000)>>11)*3;
|
|
if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
|
|
else {xx = -x - asncs.x[n]; eps=1.02; }
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
|
|
xx*asncs.x[n+7])))))+asncs.x[n+8];
|
|
t+=p;
|
|
y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
|
|
t = (m>0)?(hp1.x-t):(hp1.x+t);
|
|
res = y+t;
|
|
if (res == res+eps*((y-res)+t)) return res;
|
|
else {
|
|
r=asncs.x[n+9]+xx*asncs.x[n+10];
|
|
t=((asncs.x[n+9]-r)+xx*asncs.x[n+10])+(p+xx*asncs.x[n+11]);
|
|
if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0004; }
|
|
else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0002; }
|
|
res = p+t;
|
|
cor = (p-res)+t;
|
|
if (res == (res+eps*cor)) return res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__docos(res,z,w);
|
|
z=(w[0]-x)+w[1];
|
|
if (z>1.0e-27) return max(res,res1);
|
|
else if (z<-1.0e-27) return min(res,res1);
|
|
else return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fe80000) */
|
|
|
|
/*------------------------- 0.75 <= |x| < 0.921875 -------------*/
|
|
else
|
|
if (k < 0x3fed8000) {
|
|
n = 992+((k&0x000fe000)>>13)*13;
|
|
if (m>0) {xx = x - asncs.x[n]; eps = 1.04; }
|
|
else {xx = -x - asncs.x[n]; eps = 1.01; }
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
|
|
xx*asncs.x[n+8]))))))+asncs.x[n+9];
|
|
t+=p;
|
|
y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
|
|
t = (m>0)?(hp1.x-t):(hp1.x+t);
|
|
res = y+t;
|
|
if (res == res+eps*((y-res)+t)) return res;
|
|
else {
|
|
r=asncs.x[n+10]+xx*asncs.x[n+11];
|
|
t=((asncs.x[n+10]-r)+xx*asncs.x[n+11])+(p+xx*asncs.x[n+12]);
|
|
if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0032; }
|
|
else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0008; }
|
|
res = p+t;
|
|
cor = (p-res)+t;
|
|
if (res == (res+eps*cor)) return res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__docos(res,z,w);
|
|
z=(w[0]-x)+w[1];
|
|
if (z>1.0e-27) return max(res,res1);
|
|
else if (z<-1.0e-27) return min(res,res1);
|
|
else return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fed8000) */
|
|
|
|
/*-------------------0.921875 <= |x| < 0.953125 ------------------*/
|
|
else
|
|
if (k < 0x3fee8000) {
|
|
n = 884+((k&0x000fe000)>>13)*14;
|
|
if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
|
|
else {xx = -x - asncs.x[n]; eps =1.005; }
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
|
|
+xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
|
|
xx*asncs.x[n+9])))))))+asncs.x[n+10];
|
|
t+=p;
|
|
y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
|
|
t = (m>0)?(hp1.x-t):(hp1.x+t);
|
|
res = y+t;
|
|
if (res == res+eps*((y-res)+t)) return res;
|
|
else {
|
|
r=asncs.x[n+11]+xx*asncs.x[n+12];
|
|
t=((asncs.x[n+11]-r)+xx*asncs.x[n+12])+(p+xx*asncs.x[n+13]);
|
|
if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
|
|
else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
|
|
res = p+t;
|
|
cor = (p-res)+t;
|
|
if (res == (res+eps*cor)) return res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__docos(res,z,w);
|
|
z=(w[0]-x)+w[1];
|
|
if (z>1.0e-27) return max(res,res1);
|
|
else if (z<-1.0e-27) return min(res,res1);
|
|
else return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fee8000) */
|
|
|
|
/*--------------------0.953125 <= |x| < 0.96875 ----------------*/
|
|
else
|
|
if (k < 0x3fef0000) {
|
|
n = 768+((k&0x000fe000)>>13)*15;
|
|
if (m>0) {xx = x - asncs.x[n]; eps=1.04; }
|
|
else {xx = -x - asncs.x[n]; eps=1.005;}
|
|
t = asncs.x[n+1]*xx;
|
|
p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
|
|
xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
|
|
+xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
|
|
xx*asncs.x[n+10]))))))))+asncs.x[n+11];
|
|
t+=p;
|
|
y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
|
|
t = (m>0)?(hp1.x-t):(hp1.x+t);
|
|
res = y+t;
|
|
if (res == res+eps*((y-res)+t)) return res;
|
|
else {
|
|
r=asncs.x[n+12]+xx*asncs.x[n+13];
|
|
t=((asncs.x[n+12]-r)+xx*asncs.x[n+13])+(p+xx*asncs.x[n+14]);
|
|
if (m>0) {p = hp0.x-r; t = (((hp0.x-p)-r)-t)+hp1.x; eps=1.0030; }
|
|
else {p = hp0.x+r; t = ((hp0.x-p)+r)+(hp1.x+t); eps=1.0005; }
|
|
res = p+t;
|
|
cor = (p-res)+t;
|
|
if (res == (res+eps*cor)) return res;
|
|
else {
|
|
res1=res+1.1*cor;
|
|
z=0.5*(res1-res);
|
|
__docos(res,z,w);
|
|
z=(w[0]-x)+w[1];
|
|
if (z>1.0e-27) return max(res,res1);
|
|
else if (z<-1.0e-27) return min(res,res1);
|
|
else return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
} /* else if (k < 0x3fef0000) */
|
|
/*-----------------0.96875 <= |x| < 1 ---------------------------*/
|
|
|
|
else
|
|
if (k<0x3ff00000) {
|
|
z = 0.5*((m>0)?(1.0-x):(1.0+x));
|
|
v.x=z;
|
|
k=v.i[HIGH_HALF];
|
|
t=inroot[(k&0x001fffff)>>14]*powtwo[511-(k>>21)];
|
|
r=1.0-t*t*z;
|
|
t = t*(rt0+r*(rt1+r*(rt2+r*rt3)));
|
|
c=t*z;
|
|
t=c*(1.5-0.5*t*c);
|
|
y = (t27*c+c)-t27*c;
|
|
cc = (z-y*y)/(t+y);
|
|
p=(((((f6*z+f5)*z+f4)*z+f3)*z+f2)*z+f1)*z;
|
|
if (m<0) {
|
|
cor = (hp1.x - cc)-(y+cc)*p;
|
|
res1 = hp0.x - y;
|
|
res =res1 + cor;
|
|
if (res == res+1.002*((res1-res)+cor)) return (res+res);
|
|
else {
|
|
c=y+cc;
|
|
cc=(y-c)+cc;
|
|
__doasin(c,cc,w);
|
|
res1=hp0.x-w[0];
|
|
cor=((hp0.x-res1)-w[0])+(hp1.x-w[1]);
|
|
res = res1+cor;
|
|
cor = (res1-res)+cor;
|
|
if (res==(res+1.000001*cor)) return (res+res);
|
|
else {
|
|
res=res+res;
|
|
res1=res+1.2*cor;
|
|
return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
cor = cc+p*(y+cc);
|
|
res = y + cor;
|
|
if (res == res+1.03*((y-res)+cor)) return (res+res);
|
|
else {
|
|
c=y+cc;
|
|
cc=(y-c)+cc;
|
|
__doasin(c,cc,w);
|
|
res = w[0];
|
|
cor=w[1];
|
|
if (res==(res+1.000001*cor)) return (res+res);
|
|
else {
|
|
res=res+res;
|
|
res1=res+1.2*cor;
|
|
return __cos32(x,res,res1);
|
|
}
|
|
}
|
|
}
|
|
} /* else if (k < 0x3ff00000) */
|
|
|
|
/*---------------------------- |x|>=1 -----------------------*/
|
|
else
|
|
if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
|
|
else
|
|
if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
|
|
else {
|
|
u.i[HIGH_HALF]=0x7ff00000;
|
|
v.i[HIGH_HALF]=0x7ff00000;
|
|
u.i[LOW_HALF]=0;
|
|
v.i[LOW_HALF]=0;
|
|
return u.x/v.x;
|
|
}
|
|
}
|
|
#ifndef __ieee754_acos
|
|
strong_alias (__ieee754_acos, __acos_finite)
|
|
#endif
|