AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-10 07:10:06 +00:00
Code Issues Packages Projects Releases Wiki Activity

Format floating routines.

Browse Source
This commit is contained in:
Ondřej Bílka 2013-10-17 16:03:24 +02:00
parent e5c2c2d0c0
commit c5d5d574cb
59 changed files with 2850 additions and 2195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
ChangeLog
View File

@ -1,3 +1,64 @@
2013-10-17 Ondřej Bílka <neleai@seznam.cz>
* sysdeps/ieee754/dbl-64/dbl2mpn.c: Fix formatting.
* sysdeps/ieee754/dbl-64/dla.h: Likewise.
* sysdeps/ieee754/dbl-64/dosincos.c: Likewise.
* sysdeps/ieee754/dbl-64/e_acosh.c: Likewise.
* sysdeps/ieee754/dbl-64/e_atan2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_cosh.c: Likewise.
* sysdeps/ieee754/dbl-64/e_exp2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_fmod.c: Likewise.
* sysdeps/ieee754/dbl-64/e_gamma_r.c: Likewise.
* sysdeps/ieee754/dbl-64/e_hypot.c: Likewise.
* sysdeps/ieee754/dbl-64/e_ilogb.c: Likewise.
* sysdeps/ieee754/dbl-64/e_j0.c: Likewise.
* sysdeps/ieee754/dbl-64/e_j1.c: Likewise.
* sysdeps/ieee754/dbl-64/e_jn.c: Likewise.
* sysdeps/ieee754/dbl-64/e_log10.c: Likewise.
* sysdeps/ieee754/dbl-64/e_log2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_log.c: Likewise.
* sysdeps/ieee754/dbl-64/e_remainder.c: Likewise.
* sysdeps/ieee754/dbl-64/e_rem_pio2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_sinh.c: Likewise.
* sysdeps/ieee754/dbl-64/e_sqrt.c: Likewise.
* sysdeps/ieee754/dbl-64/halfulp.c: Likewise.
* sysdeps/ieee754/dbl-64/k_rem_pio2.c: Likewise.
* sysdeps/ieee754/dbl-64/MathLib.h: Likewise.
* sysdeps/ieee754/dbl-64/mpa-arch.h: Likewise.
* sysdeps/ieee754/dbl-64/mpa.c: Likewise.
* sysdeps/ieee754/dbl-64/mpatan.c: Likewise.
* sysdeps/ieee754/dbl-64/mpn2dbl.c: Likewise.
* sysdeps/ieee754/dbl-64/mptan.c: Likewise.
* sysdeps/ieee754/dbl-64/mydefs.h: Likewise.
* sysdeps/ieee754/dbl-64/s_asinh.c: Likewise.
* sysdeps/ieee754/dbl-64/s_atan.c: Likewise.
* sysdeps/ieee754/dbl-64/s_cbrt.c: Likewise.
* sysdeps/ieee754/dbl-64/s_ceil.c: Likewise.
* sysdeps/ieee754/dbl-64/s_copysign.c: Likewise.
* sysdeps/ieee754/dbl-64/s_erf.c: Likewise.
* sysdeps/ieee754/dbl-64/s_expm1.c: Likewise.
* sysdeps/ieee754/dbl-64/s_fabs.c: Likewise.
* sysdeps/ieee754/dbl-64/s_finite.c: Likewise.
* sysdeps/ieee754/dbl-64/s_floor.c: Likewise.
* sysdeps/ieee754/dbl-64/s_frexp.c: Likewise.
* sysdeps/ieee754/dbl-64/s_isinf.c: Likewise.
* sysdeps/ieee754/dbl-64/s_isinf_ns.c: Likewise.
* sysdeps/ieee754/dbl-64/s_isnan.c: Likewise.
* sysdeps/ieee754/dbl-64/s_llround.c: Likewise.
* sysdeps/ieee754/dbl-64/s_log1p.c: Likewise.
* sysdeps/ieee754/dbl-64/s_logb.c: Likewise.
* sysdeps/ieee754/dbl-64/s_lrint.c: Likewise.
* sysdeps/ieee754/dbl-64/s_modf.c: Likewise.
* sysdeps/ieee754/dbl-64/s_nearbyint.c: Likewise.
* sysdeps/ieee754/dbl-64/s_remquo.c: Likewise.
* sysdeps/ieee754/dbl-64/s_rint.c: Likewise.
* sysdeps/ieee754/dbl-64/s_scalbln.c: Likewise.
* sysdeps/ieee754/dbl-64/s_scalbn.c: Likewise.
* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c: Likewise.
* sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
* sysdeps/ieee754/dbl-64/s_tanh.c: Likewise.
2013-10-17 Joseph Myers <joseph@codesourcery.com>
[BZ #16041]

30
sysdeps/ieee754/dbl-64/MathLib.h
View File

@ -33,14 +33,14 @@
/* It returns the original status of these modes. */
/* See further explanations of usage in DPChange.h */
/********************************************************************/
unsigned short Init_Lib(void);
unsigned short Init_Lib (void);
/********************************************************************/
/* Function that changes the precision and rounding modes to the */
/* specified by the argument received. See further explanations in */
/* DPChange.h */
/********************************************************************/
void Exit_Lib(unsigned short);
void Exit_Lib (unsigned short);
/* The asin() function calculates the arc sine of its argument. */
@ -48,7 +48,7 @@ void Exit_Lib(unsigned short);
/* (between -PI/2 and PI/2). */
/* If the argument is greater than 1 or less than -1 it returns */
/* a NaN. */
double uasin(double );
double uasin (double);
/* The acos() function calculates the arc cosine of its argument. */
@ -56,47 +56,45 @@ double uasin(double );
/* (between -PI/2 and PI/2). */
/* If the argument is greater than 1 or less than -1 it returns */
/* a NaN. */
double uacos(double );
double uacos (double);
/* The atan() function calculates the arctanget of its argument. */
/* The function returns the arc tangent in radians */
/* (between -PI/2 and PI/2). */
double uatan(double );
double uatan (double);
/* The uatan2() function calculates the arc tangent of the two arguments x */
/* and y (x is the right argument and y is the left one).The signs of both */
/* arguments are used to determine the quadrant of the result. */
/* The function returns the result in radians, which is between -PI and PI */
double uatan2(double ,double );
double uatan2 (double, double);
/* Compute log(x). The base of log is e (natural logarithm) */
double ulog(double );
double ulog (double);
/* Compute e raised to the power of argument x. */
double uexp(double );
double uexp (double);
/* Compute sin(x). The argument x is assumed to be given in radians.*/
double usin(double );
double usin (double);
/* Compute cos(x). The argument x is assumed to be given in radians.*/
double ucos(double );
double ucos (double);
/* Compute tan(x). The argument x is assumed to be given in radians.*/
double utan(double );
double utan (double);
/* Compute the square root of non-negative argument x. */
/* If x is negative the returned value is NaN. */
double usqrt(double );
double usqrt (double);
/* Compute x raised to the power of y, where x is the left argument */
/* and y is the right argument. The function returns a NaN if x<0. */
/* If x equals zero it returns -inf */
double upow(double , double );
double upow (double, double);
/* Computing x mod y, where x is the left argument and y is the */
/* right one. */
double uremainder(double , double );
double uremainder (double, double);
#endif

8
sysdeps/ieee754/dbl-64/dbl2mpn.c
View File

@ -40,14 +40,14 @@ __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
#if BITS_PER_MP_LIMB == 32
res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction. */
res_ptr[1] = u.ieee.mantissa0; /* High-order 20 bits. */
#define N 2
# define N 2
#elif BITS_PER_MP_LIMB == 64
/* Hopefully the compiler will combine the two bitfield extracts
and this composition into just the original quadword extract. */
res_ptr[0] = ((mp_limb_t) u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
#define N 1
# define N 1
#else
#error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
# error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
#endif
/* The format does not fill the last limb. There are some zeros. */
#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
@ -73,7 +73,7 @@ __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
#if N == 2
res_ptr[N - 1] = res_ptr[1] << cnt
| (N - 1)
* (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
* (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
res_ptr[0] <<= cnt;
#else
res_ptr[N - 1] <<= cnt;

86
sysdeps/ieee754/dbl-64/dla.h
View File

@ -61,13 +61,13 @@
/* storage variables of type double. */
#ifdef DLA_FMS
# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \
z=x*y; zz=DLA_FMS(x,y,z);
# define EMULV(x, y, z, zz, p, hx, tx, hy, ty) \
z = x * y; zz = DLA_FMS (x, y, z);
#else
# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \
p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \
p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \
z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty;
# define EMULV(x, y, z, zz, p, hx, tx, hy, ty) \
p = CN * (x); hx = ((x) - p) + p; tx = (x) - hx; \
p = CN * (y); hy = ((y) - p) + p; ty = (y) - hy; \
z = (x) * (y); zz = (((hx * hy - z) + hx * ty) + tx * hy) + tx * ty;
#endif
@ -93,11 +93,11 @@
/* are assumed to be double-length numbers. r,s are temporary */
/* storage variables of type double. */
#define ADD2(x,xx,y,yy,z,zz,r,s) \
r=(x)+(y); s=(ABS(x)>ABS(y)) ? \
(((((x)-r)+(y))+(yy))+(xx)) : \
(((((y)-r)+(x))+(xx))+(yy)); \
z=r+s; zz=(r-z)+s;
#define ADD2(x, xx, y, yy, z, zz, r, s) \
r = (x) + (y); s = (ABS (x) > ABS (y)) ? \
(((((x) - r) + (y)) + (yy)) + (xx)) : \
(((((y) - r) + (x)) + (xx)) + (yy)); \
z = r + s; zz = (r - z) + s;
/* Double-length subtraction, Dekker. The macro produces a double-length */
@ -106,11 +106,11 @@
/* are assumed to be double-length numbers. r,s are temporary */
/* storage variables of type double. */
#define SUB2(x,xx,y,yy,z,zz,r,s) \
r=(x)-(y); s=(ABS(x)>ABS(y)) ? \
(((((x)-r)-(y))-(yy))+(xx)) : \
((((x)-((y)+r))+(xx))-(yy)); \
z=r+s; zz=(r-z)+s;
#define SUB2(x, xx, y, yy, z, zz, r, s) \
r = (x) - (y); s = (ABS (x) > ABS (y)) ? \
(((((x) - r) - (y)) - (yy)) + (xx)) : \
((((x) - ((y) + r)) + (xx)) - (yy)); \
z = r + s; zz = (r - z) + s;
/* Double-length multiplication, Dekker. The macro produces a double-length */
@ -119,9 +119,9 @@
/* are assumed to be double-length numbers. p,hx,tx,hy,ty,q,c,cc are */
/* temporary storage variables of type double. */
#define MUL2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc) \
MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \
cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc;
#define MUL2(x, xx, y, yy, z, zz, p, hx, tx, hy, ty, q, c, cc) \
MUL12 (x, y, c, cc, p, hx, tx, hy, ty, q) \
cc = ((x) * (yy) + (xx) * (y)) + cc; z = c + cc; zz = (c - z) + cc;
/* Double-length division, Dekker. The macro produces a double-length */
@ -142,18 +142,18 @@
/* are assumed to be double-length numbers. r,rr,s,ss,u,uu,w */
/* are temporary storage variables of type double. */
#define ADD2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \
r=(x)+(y); \
if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \
else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \
if (rr!=0.0) { \
z=r+s; zz=(r-z)+s; } \
else { \
ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \
u=r+s; \
uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \
w=uu+ss; z=u+w; \
zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; }
#define ADD2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \
r = (x) + (y); \
if (ABS (x) > ABS (y)) { rr = ((x) - r) + (y); s = (rr + (yy)) + (xx); } \
else { rr = ((y) - r) + (x); s = (rr + (xx)) + (yy); } \
if (rr != 0.0) { \
z = r + s; zz = (r - z) + s; } \
else { \
ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) + (yy)) : (((yy) - s) + (xx));\
u = r + s; \
uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \
w = uu + ss; z = u + w; \
zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); }
/* Double-length subtraction, slower but more accurate than SUB2. */
@ -163,15 +163,15 @@
/* are assumed to be double-length numbers. r,rr,s,ss,u,uu,w */
/* are temporary storage variables of type double. */
#define SUB2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \
r=(x)-(y); \
if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \
else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \
if (rr!=0.0) { \
z=r+s; zz=(r-z)+s; } \
else { \
ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \
u=r+s; \
uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \
w=uu+ss; z=u+w; \
zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; }
#define SUB2A(x, xx, y, yy, z, zz, r, rr, s, ss, u, uu, w) \
r = (x) - (y); \
if (ABS (x) > ABS (y)) { rr = ((x) - r) - (y); s = (rr - (yy)) + (xx); } \
else { rr = (x) - ((y) + r); s = (rr + (xx)) - (yy); } \
if (rr != 0.0) { \
z = r + s; zz = (r - z) + s; } \
else { \
ss = (ABS (xx) > ABS (yy)) ? (((xx) - s) - (yy)) : ((xx) - ((yy) + s)); \
u = r + s; \
uu = (ABS (r) > ABS (s)) ? ((r - u) + s) : ((s - u) + r); \
w = uu + ss; z = u + w; \
zz = (ABS (u) > ABS (w)) ? ((u - z) + w) : ((w - z) + u); }

240
sysdeps/ieee754/dbl-64/dosincos.c
View File

@ -57,49 +57,52 @@ extern const union
void
SECTION
__dubsin(double x, double dx, double v[]) {
double r,s,c,cc,d,dd,d2,dd2,e,ee,
sn,ssn,cs,ccs,ds,dss,dc,dcc;
__dubsin (double x, double dx, double v[])
{
double r, s, c, cc, d, dd, d2, dd2, e, ee,
sn, ssn, cs, ccs, ds, dss, dc, dcc;
#ifndef DLA_FMS
double p,hx,tx,hy,ty,q;
double p, hx, tx, hy, ty, q;
#endif
mynumber u;
int4 k;
u.x=x+big.x;
k = u.i[LOW_HALF]<<2;
x=x-(u.x-big.x);
d=x+dx;
dd=(x-d)+dx;
/* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
sn=__sincostab.x[k]; /* */
ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */
cs=__sincostab.x[k+2]; /* */
ccs=__sincostab.x[k+3]; /* */
MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */
ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */
ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* for sin */
MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,d,dd,ds,dss,r,s); /* ds=sin(t) */
u.x = x + big.x;
k = u.i[LOW_HALF] << 2;
x = x - (u.x - big.x);
d = x + dx;
dd = (x - d) + dx;
/* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
sn = __sincostab.x[k]; /* */
ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */
cs = __sincostab.x[k + 2]; /* */
ccs = __sincostab.x[k + 3]; /* */
/* Taylor series for sin ds=sin(t) */
MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, d, dd, ds, dss, r, s);
MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc); ;/* Taylor */
ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* series */
ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* for cos */
ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc); /* dc=cos(t) */
/* Taylor series for cos dc=cos(t) */
MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
SUB2(e,ee,dc,dcc,e,ee,r,s);
ADD2(e,ee,sn,ssn,e,ee,r,s); /* e+ee=sin(x+dx) */
MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
SUB2 (e, ee, dc, dcc, e, ee, r, s);
ADD2 (e, ee, sn, ssn, e, ee, r, s); /* e+ee=sin(x+dx) */
v[0]=e;
v[1]=ee;
v[0] = e;
v[1] = ee;
}
/**********************************************************************/
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
@ -110,64 +113,65 @@ __dubsin(double x, double dx, double v[]) {
void
SECTION
__dubcos(double x, double dx, double v[]) {
double r,s,c,cc,d,dd,d2,dd2,e,ee,
sn,ssn,cs,ccs,ds,dss,dc,dcc;
__dubcos (double x, double dx, double v[])
{
double r, s, c, cc, d, dd, d2, dd2, e, ee,
sn, ssn, cs, ccs, ds, dss, dc, dcc;
#ifndef DLA_FMS
double p,hx,tx,hy,ty,q;
double p, hx, tx, hy, ty, q;
#endif
mynumber u;
int4 k;
u.x=x+big.x;
k = u.i[LOW_HALF]<<2;
x=x-(u.x-big.x);
d=x+dx;
dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
sn=__sincostab.x[k]; /* */
ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */
cs=__sincostab.x[k+2]; /* */
ccs=__sincostab.x[k+3]; /* */
MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,d,dd,ds,dss,r,s);
u.x = x + big.x;
k = u.i[LOW_HALF] << 2;
x = x - (u.x - big.x);
d = x + dx;
dd = (x - d) + dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */
MUL2 (d, dd, d, dd, d2, dd2, p, hx, tx, hy, ty, q, c, cc);
sn = __sincostab.x[k]; /* */
ssn = __sincostab.x[k + 1]; /* sin(Xi) and cos(Xi) */
cs = __sincostab.x[k + 2]; /* */
ccs = __sincostab.x[k + 3]; /* */
MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, d, dd, ds, dss, r, s);
MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
MUL2(cs,ccs,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
MUL2(dc,dcc,sn,ssn,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
MUL2 (cs, ccs, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
MUL2 (dc, dcc, sn, ssn, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s);
MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,s3.x,ss3.x,ds,dss,r,s);
MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
MUL2(d,dd,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc);
ADD2(ds,dss,d,dd,ds,dss,r,s);
MUL2(d2,dd2,c8.x,cc8.x,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(dc,dcc,c6.x,cc6.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(dc,dcc,c4.x,cc4.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(dc,dcc,c2.x,cc2.x,dc,dcc,r,s);
MUL2(d2,dd2,dc,dcc,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
MUL2(sn,ssn,ds,dss,e,ee,p,hx,tx,hy,ty,q,c,cc);
MUL2(dc,dcc,cs,ccs,dc,dcc,p,hx,tx,hy,ty,q,c,cc);
ADD2(e,ee,dc,dcc,e,ee,r,s);
SUB2(cs,ccs,e,ee,e,ee,r,s);
MUL2 (d2, dd2, s7.x, ss7.x, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, s5.x, ss5.x, ds, dss, r, s);
MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, s3.x, ss3.x, ds, dss, r, s);
MUL2 (d2, dd2, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
MUL2 (d, dd, ds, dss, ds, dss, p, hx, tx, hy, ty, q, c, cc);
ADD2 (ds, dss, d, dd, ds, dss, r, s);
MUL2 (d2, dd2, c8.x, cc8.x, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c6.x, cc6.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c4.x, cc4.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (dc, dcc, c2.x, cc2.x, dc, dcc, r, s);
MUL2 (d2, dd2, dc, dcc, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
MUL2 (sn, ssn, ds, dss, e, ee, p, hx, tx, hy, ty, q, c, cc);
MUL2 (dc, dcc, cs, ccs, dc, dcc, p, hx, tx, hy, ty, q, c, cc);
ADD2 (e, ee, dc, dcc, e, ee, r, s);
SUB2 (cs, ccs, e, ee, e, ee, r, s);
v[0]=e;
v[1]=ee;
v[0] = e;
v[1] = ee;
}
/**********************************************************************/
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
@ -175,29 +179,45 @@ __dubcos(double x, double dx, double v[]) {
/**********************************************************************/
void
SECTION
__docos(double x, double dx, double v[]) {
double y,yy,p,w[2];
if (x>0) {y=x; yy=dx;}
else {y=-x; yy=-dx;}
if (y<0.5*hp0.x) /* y< PI/4 */
{__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
else if (y<1.5*hp0.x) { /* y< 3/4 * PI */
p=hp0.x-y; /* p = PI/2 - y */
yy=hp1.x-yy;
y=p+yy;
yy=(p-y)+yy;
if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
/* cos(x) = sin ( 90 - x ) */
else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
}
}
else { /* y>= 3/4 * PI */
p=2.0*hp0.x-y; /* p = PI- y */
yy=2.0*hp1.x-yy;
y=p+yy;
yy=(p-y)+yy;
__dubcos(y,yy,w);
v[0]=-w[0];
v[1]=-w[1];
}
__docos (double x, double dx, double v[])
{
double y, yy, p, w[2];
if (x > 0)
{
y = x; yy = dx;
}
else
{
y = -x; yy = -dx;
}
if (y < 0.5 * hp0.x) /* y< PI/4 */
{
__dubcos (y, yy, w); v[0] = w[0]; v[1] = w[1];
}
else if (y < 1.5 * hp0.x) /* y< 3/4 * PI */
{
p = hp0.x - y; /* p = PI/2 - y */
yy = hp1.x - yy;
y = p + yy;
yy = (p - y) + yy;
if (y > 0)
{
__dubsin (y, yy, w); v[0] = w[0]; v[1] = w[1];
}
/* cos(x) = sin ( 90 - x ) */
else
{
__dubsin (-y, -yy, w); v[0] = -w[0]; v[1] = -w[1];
}
}
else /* y>= 3/4 * PI */
{
p = 2.0 * hp0.x - y; /* p = PI- y */
yy = 2.0 * hp1.x - yy;
y = p + yy;
yy = (p - y) + yy;
__dubcos (y, yy, w);
v[0] = -w[0];
v[1] = -w[1];
}
}

55
sysdeps/ieee754/dbl-64/e_acosh.c
View File

@ -28,31 +28,42 @@
#include <math_private.h>
static const double
one = 1.0,
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
one = 1.0,
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double
__ieee754_acosh(double x)
__ieee754_acosh (double x)
{
double t;
int32_t hx;
u_int32_t lx;
EXTRACT_WORDS(hx,lx,x);
if(hx<0x3ff00000) { /* x < 1 */
return (x-x)/(x-x);
} else if(hx >=0x41b00000) { /* x > 2**28 */
if(hx >=0x7ff00000) { /* x is inf of NaN */
return x+x;
} else
return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */
} else if(((hx-0x3ff00000)|lx)==0) {
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t=x*x;
return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
} else { /* 1<x<2 */
t = x-one;
return __log1p(t+__ieee754_sqrt(2.0*t+t*t));
double t;
int32_t hx;
u_int32_t lx;
EXTRACT_WORDS (hx, lx, x);
if (hx < 0x3ff00000) /* x < 1 */
{
return (x - x) / (x - x);
}
else if (hx >= 0x41b00000) /* x > 2**28 */
{
if (hx >= 0x7ff00000) /* x is inf of NaN */
{
return x + x;
}
else
return __ieee754_log (x) + ln2; /* acosh(huge)=log(2x) */
}
else if (((hx - 0x3ff00000) | lx) == 0)
{
return 0.0; /* acosh(1) = 0 */
}
else if (hx > 0x40000000) /* 2**28 > x > 2 */
{
t = x * x;
return __ieee754_log (2.0 * x - one / (x + __ieee754_sqrt (t - one)));
}
else /* 1<x<2 */
{
t = x - one;
return __log1p (t + __ieee754_sqrt (2.0 * t + t * t));
}
}
strong_alias (__ieee754_acosh, __acosh_finite)

20
sysdeps/ieee754/dbl-64/e_atan2.c
View File

@ -72,14 +72,14 @@ __ieee754_atan2 (double y, double x)
int i, de, ux, dx, uy, dy;
static const int pr[MM] = { 6, 8, 10, 20, 32 };
double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8,
z, zz, cor, s1, ss1, s2, ss2;
z, zz, cor, s1, ss1, s2, ss2;
#ifndef DLA_FMS
double t4, t5, t6;
#endif
number num;
static const int ep = 59768832, /* 57*16**5 */
em = -59768832; /* -57*16**5 */
static const int ep = 59768832, /* 57*16**5 */
em = -59768832; /* -57*16**5 */
/* x=NaN or y=NaN */
num.d = x;
@ -294,17 +294,17 @@ __ieee754_atan2 (double y, double x)
if (i < 112)
{
if (i < 48)
u9 = u91.d; /* u < 1/4 */
u9 = u91.d; /* u < 1/4 */
else
u9 = u92.d;
} /* 1/4 <= u < 1/2 */
} /* 1/4 <= u < 1/2 */
else
{
if (i < 176)
u9 = u93.d; /* 1/2 <= u < 3/4 */
u9 = u93.d; /* 1/2 <= u < 3/4 */
else
u9 = u94.d;
} /* 3/4 <= u <= 1 */
} /* 3/4 <= u <= 1 */
if ((z = t1 + (zz - u9 * t1)) == t1 + (zz + u9 * t1))
return signArctan2 (y, z);
@ -312,9 +312,9 @@ __ieee754_atan2 (double y, double x)
EADD (t1, du, v, vv);
s1 = v * (hij[i][11].d
+ v * (hij[i][12].d
+ v * (hij[i][13].d
+ v * (hij[i][14].d
+ v * hij[i][15].d))));
+ v * (hij[i][13].d
+ v * (hij[i][14].d
+ v * hij[i][15].d))));
ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);

72
sysdeps/ieee754/dbl-64/e_cosh.c
View File

@ -34,49 +34,55 @@
#include <math.h>
#include <math_private.h>
static const double one = 1.0, half=0.5, huge = 1.0e300;
static const double one = 1.0, half = 0.5, huge = 1.0e300;
double
__ieee754_cosh (double x)
{
double t,w;
int32_t ix;
u_int32_t lx;
double t, w;
int32_t ix;
u_int32_t lx;
/* High word of |x|. */
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
/* High word of |x|. */
GET_HIGH_WORD (ix, x);
ix &= 0x7fffffff;
/* |x| in [0,22] */
if (ix < 0x40360000) {
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if(ix<0x3fd62e43) {
t = __expm1(fabs(x));
w = one+t;
if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */
return one+(t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
t = __ieee754_exp(fabs(x));
return half*t+half/t;
/* |x| in [0,22] */
if (ix < 0x40360000)
{
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43)
{
t = __expm1 (fabs (x));
w = one + t;
if (ix < 0x3c800000)
return w; /* cosh(tiny) = 1 */
return one + (t * t) / (w + w);
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862e42) return half*__ieee754_exp(fabs(x));
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
t = __ieee754_exp (fabs (x));
return half * t + half / t;
}
/* |x| in [log(maxdouble), overflowthresold] */
GET_LOW_WORD(lx,x);
if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
w = __ieee754_exp(half*fabs(x));
t = half*w;
return t*w;
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862e42)
return half * __ieee754_exp (fabs (x));
/* x is INF or NaN */
if(ix>=0x7ff00000) return x*x;
/* |x| in [log(maxdouble), overflowthresold] */
GET_LOW_WORD (lx, x);
if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (u_int32_t) 0x8fb9f87d)))
{
w = __ieee754_exp (half * fabs (x));
t = half * w;
return t * w;
}
/* |x| > overflowthresold, cosh(x) overflow */
return huge*huge;
/* x is INF or NaN */
if (ix >= 0x7ff00000)
return x * x;
/* |x| > overflowthresold, cosh(x) overflow */
return huge * huge;
}
strong_alias (__ieee754_cosh, __cosh_finite)

8
sysdeps/ieee754/dbl-64/e_exp2.c
View File

@ -46,7 +46,7 @@ __ieee754_exp2 (double x)
if (__builtin_expect (isless (x, himark), 1))
{
/* Exceptional cases: */
if (__builtin_expect (! isgreaterequal (x, lomark), 0))
if (__builtin_expect (!isgreaterequal (x, lomark), 0))
{
if (__isinf (x))
/* e^-inf == 0, with no error. */
@ -93,7 +93,7 @@ __ieee754_exp2 (double x)
/* 3. Compute ex2 = 2^(t/512+e+ex). */
ex2_u.d = exp2_accuratetable[tval & 511];
tval >>= 9;
unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
unsafe = abs (tval) >= -DBL_MIN_EXP - 1;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.d = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
@ -106,7 +106,7 @@ __ieee754_exp2 (double x)
* x + .055504110254308625)
* x + .240226506959100583)
* x + .69314718055994495) * ex2_u.d;
math_opt_barrier (x22);
math_opt_barrier (x22);
}
/* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */
@ -119,6 +119,6 @@ __ieee754_exp2 (double x)
}
else
/* Return x, if x is a NaN or Inf; or overflow, otherwise. */
return TWO1023*x;
return TWO1023 * x;
}
strong_alias (__ieee754_exp2, __exp2_finite)

233
sysdeps/ieee754/dbl-64/e_fmod.c
View File

@ -18,111 +18,156 @@
#include <math.h>
#include <math_private.h>
static const double one = 1.0, Zero[] = {0.0, -0.0,};
static const double one = 1.0, Zero[] = { 0.0, -0.0, };
double
__ieee754_fmod (double x, double y)
{
int32_t n,hx,hy,hz,ix,iy,sx,i;
u_int32_t lx,ly,lz;
int32_t n, hx, hy, hz, ix, iy, sx, i;
u_int32_t lx, ly, lz;
EXTRACT_WORDS(hx,lx,x);
EXTRACT_WORDS(hy,ly,y);
sx = hx&0x80000000; /* sign of x */
hx ^=sx; /* |x| */
hy &= 0x7fffffff; /* |y| */
EXTRACT_WORDS (hx, lx, x);
EXTRACT_WORDS (hy, ly, y);
sx = hx & 0x80000000; /* sign of x */
hx ^= sx; /* |x| */
hy &= 0x7fffffff; /* |y| */
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
return (x*y)/(x*y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
if(lx==ly)
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
/* purge off exception values */
if ((hy | ly) == 0 || (hx >= 0x7ff00000) || /* y=0,or x not finite */
((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */
return (x * y) / (x * y);
if (hx <= hy)
{
if ((hx < hy) || (lx < ly))
return x; /* |x|<|y| return x */
if (lx == ly)
return Zero[(u_int32_t) sx >> 31]; /* |x|=|y| return x*0*/
}
/* determine ix = ilogb(x) */
if (__builtin_expect (hx < 0x00100000, 0)) /* subnormal x */
{
if (hx == 0)
{
for (ix = -1043, i = lx; i > 0; i <<= 1)
ix -= 1;
}
/* determine ix = ilogb(x) */
if(__builtin_expect(hx<0x00100000, 0)) { /* subnormal x */
if(hx==0) {
for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
} else {
for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
}
} else ix = (hx>>20)-1023;
/* determine iy = ilogb(y) */
if(__builtin_expect(hy<0x00100000, 0)) { /* subnormal y */
if(hy==0) {
for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
} else {
for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
}
} else iy = (hy>>20)-1023;
/* set up {hx,lx}, {hy,ly} and align y to x */
if(__builtin_expect(ix >= -1022, 1))
hx = 0x00100000|(0x000fffff&hx);
else { /* subnormal x, shift x to normal */
n = -1022-ix;
if(n<=31) {
hx = (hx<<n)|(lx>>(32-n));
lx <<= n;
} else {
hx = lx<<(n-32);
lx = 0;
}
}
if(__builtin_expect(iy >= -1022, 1))
hy = 0x00100000|(0x000fffff&hy);
else { /* subnormal y, shift y to normal */
n = -1022-iy;
if(n<=31) {
hy = (hy<<n)|(ly>>(32-n));
ly <<= n;
} else {
hy = ly<<(n-32);
ly = 0;
}
else
{
for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
ix -= 1;
}
}
else
ix = (hx >> 20) - 1023;
/* fix point fmod */
n = ix - iy;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
else {
if((hz|lz)==0) /* return sign(x)*0 */
return Zero[(u_int32_t)sx>>31];
hx = hz+hz+(lz>>31); lx = lz+lz;
}
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) /* return sign(x)*0 */
return Zero[(u_int32_t)sx>>31];
while(hx<0x00100000) { /* normalize x */
hx = hx+hx+(lx>>31); lx = lx+lx;
/* determine iy = ilogb(y) */
if (__builtin_expect (hy < 0x00100000, 0)) /* subnormal y */
{
if (hy == 0)
{
for (iy = -1043, i = ly; i > 0; i <<= 1)
iy -= 1;
}
if(__builtin_expect(iy>= -1022, 1)) { /* normalize output */
hx = ((hx-0x00100000)|((iy+1023)<<20));
INSERT_WORDS(x,hx|sx,lx);
} else { /* subnormal output */
n = -1022 - iy;
if(n<=20) {
lx = (lx>>n)|((u_int32_t)hx<<(32-n));
hx >>= n;
} else if (n<=31) {
lx = (hx<<(32-n))|(lx>>n); hx = sx;
} else {
lx = hx>>(n-32); hx = sx;
}
INSERT_WORDS(x,hx|sx,lx);
x *= one; /* create necessary signal */
else
{
for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
iy -= 1;
}
return x; /* exact output */
}
else
iy = (hy >> 20) - 1023;
/* set up {hx,lx}, {hy,ly} and align y to x */
if (__builtin_expect (ix >= -1022, 1))
hx = 0x00100000 | (0x000fffff & hx);
else /* subnormal x, shift x to normal */
{
n = -1022 - ix;
if (n <= 31)
{
hx = (hx << n) | (lx >> (32 - n));
lx <<= n;
}
else
{
hx = lx << (n - 32);
lx = 0;
}
}
if (__builtin_expect (iy >= -1022, 1))
hy = 0x00100000 | (0x000fffff & hy);
else /* subnormal y, shift y to normal */
{
n = -1022 - iy;
if (n <= 31)
{
hy = (hy << n) | (ly >> (32 - n));
ly <<= n;
}
else
{
hy = ly << (n - 32);
ly = 0;
}
}
/* fix point fmod */
n = ix - iy;
while (n--)
{
hz = hx - hy; lz = lx - ly; if (lx < ly)
hz -= 1;
if (hz < 0)
{
hx = hx + hx + (lx >> 31); lx = lx + lx;
}
else
{
if ((hz | lz) == 0) /* return sign(x)*0 */
return Zero[(u_int32_t) sx >> 31];
hx = hz + hz + (lz >> 31); lx = lz + lz;
}
}
hz = hx - hy; lz = lx - ly; if (lx < ly)
hz -= 1;
if (hz >= 0)
{
hx = hz; lx = lz;
}
/* convert back to floating value and restore the sign */
if ((hx | lx) == 0) /* return sign(x)*0 */
return Zero[(u_int32_t) sx >> 31];
while (hx < 0x00100000) /* normalize x */
{
hx = hx + hx + (lx >> 31); lx = lx + lx;
iy -= 1;
}
if (__builtin_expect (iy >= -1022, 1)) /* normalize output */
{
hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
INSERT_WORDS (x, hx | sx, lx);
}
else /* subnormal output */
{
n = -1022 - iy;
if (n <= 20)
{
lx = (lx >> n) | ((u_int32_t) hx << (32 - n));
hx >>= n;
}
else if (n <= 31)
{
lx = (hx << (32 - n)) | (lx >> n); hx = sx;
}
else
{
lx = hx >> (n - 32); hx = sx;
}
INSERT_WORDS (x, hx | sx, lx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
}
strong_alias (__ieee754_fmod, __fmod_finite)

2
sysdeps/ieee754/dbl-64/e_gamma_r.c
View File

@ -135,7 +135,7 @@ __ieee754_gamma_r (double x, int *signgamp)
*signgamp = 0;
return (x - x) / (x - x);
}
if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx==0, 0))
if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx == 0, 0))
{
/* x == -Inf. According to ISO this is NaN. */
*signgamp = 0;

157
sysdeps/ieee754/dbl-64/e_hypot.c
View File

@ -46,76 +46,101 @@
#include <math_private.h>
double
__ieee754_hypot(double x, double y)
__ieee754_hypot (double x, double y)
{
double a,b,t1,t2,y1,y2,w;
int32_t j,k,ha,hb;
double a, b, t1, t2, y1, y2, w;
int32_t j, k, ha, hb;
GET_HIGH_WORD(ha,x);
ha &= 0x7fffffff;
GET_HIGH_WORD(hb,y);
hb &= 0x7fffffff;
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
SET_HIGH_WORD(a,ha); /* a <- |a| */
SET_HIGH_WORD(b,hb); /* b <- |b| */
if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
k=0;
if(__builtin_expect(ha > 0x5f300000, 0)) { /* a>2**500 */
if(ha >= 0x7ff00000) { /* Inf or NaN */
u_int32_t low;
w = a+b; /* for sNaN */
GET_LOW_WORD(low,a);
if(((ha&0xfffff)|low)==0) w = a;
GET_LOW_WORD(low,b);
if(((hb^0x7ff00000)|low)==0) w = b;
return w;
}
/* scale a and b by 2**-600 */
ha -= 0x25800000; hb -= 0x25800000; k += 600;
SET_HIGH_WORD(a,ha);
SET_HIGH_WORD(b,hb);
GET_HIGH_WORD (ha, x);
ha &= 0x7fffffff;
GET_HIGH_WORD (hb, y);
hb &= 0x7fffffff;
if (hb > ha)
{
a = y; b = x; j = ha; ha = hb; hb = j;
}
else
{
a = x; b = y;
}
SET_HIGH_WORD (a, ha); /* a <- |a| */
SET_HIGH_WORD (b, hb); /* b <- |b| */
if ((ha - hb) > 0x3c00000)
{
return a + b;
} /* x/y > 2**60 */
k = 0;
if (__builtin_expect (ha > 0x5f300000, 0)) /* a>2**500 */
{
if (ha >= 0x7ff00000) /* Inf or NaN */
{
u_int32_t low;
w = a + b; /* for sNaN */
GET_LOW_WORD (low, a);
if (((ha & 0xfffff) | low) == 0)
w = a;
GET_LOW_WORD (low, b);
if (((hb ^ 0x7ff00000) | low) == 0)
w = b;
return w;
}
if(__builtin_expect(hb < 0x20b00000, 0)) { /* b < 2**-500 */
if(hb <= 0x000fffff) { /* subnormal b or 0 */
u_int32_t low;
GET_LOW_WORD(low,b);
if((hb|low)==0) return a;
t1=0;
SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */
b *= t1;
a *= t1;
k -= 1022;
} else { /* scale a and b by 2^600 */
ha += 0x25800000; /* a *= 2^600 */
hb += 0x25800000; /* b *= 2^600 */
k -= 600;
SET_HIGH_WORD(a,ha);
SET_HIGH_WORD(b,hb);
}
/* scale a and b by 2**-600 */
ha -= 0x25800000; hb -= 0x25800000; k += 600;
SET_HIGH_WORD (a, ha);
SET_HIGH_WORD (b, hb);
}
if (__builtin_expect (hb < 0x20b00000, 0)) /* b < 2**-500 */
{
if (hb <= 0x000fffff) /* subnormal b or 0 */
{
u_int32_t low;
GET_LOW_WORD (low, b);
if ((hb | low) == 0)
return a;
t1 = 0;
SET_HIGH_WORD (t1, 0x7fd00000); /* t1=2^1022 */
b *= t1;
a *= t1;
k -= 1022;
}
/* medium size a and b */
w = a-b;
if (w>b) {
t1 = 0;
SET_HIGH_WORD(t1,ha);
t2 = a-t1;
w = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
y1 = 0;
SET_HIGH_WORD(y1,hb);
y2 = b - y1;
t1 = 0;
SET_HIGH_WORD(t1,ha+0x00100000);
t2 = a - t1;
w = __ieee754_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
else /* scale a and b by 2^600 */
{
ha += 0x25800000; /* a *= 2^600 */
hb += 0x25800000; /* b *= 2^600 */
k -= 600;
SET_HIGH_WORD (a, ha);
SET_HIGH_WORD (b, hb);
}
if(k!=0) {
u_int32_t high;
t1 = 1.0;
GET_HIGH_WORD(high,t1);
SET_HIGH_WORD(t1,high+(k<<20));
return t1*w;
} else return w;
}
/* medium size a and b */
w = a - b;
if (w > b)
{
t1 = 0;
SET_HIGH_WORD (t1, ha);
t2 = a - t1;
w = __ieee754_sqrt (t1 * t1 - (b * (-b) - t2 * (a + t1)));
}
else
{
a = a + a;
y1 = 0;
SET_HIGH_WORD (y1, hb);
y2 = b - y1;
t1 = 0;
SET_HIGH_WORD (t1, ha + 0x00100000);
t2 = a - t1;
w = __ieee754_sqrt (t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)));
}
if (k != 0)
{
u_int32_t high;
t1 = 1.0;
GET_HIGH_WORD (high, t1);
SET_HIGH_WORD (t1, high + (k << 20));
return t1 * w;
}
else
return w;
}
strong_alias (__ieee754_hypot, __hypot_finite)

53
sysdeps/ieee754/dbl-64/e_ilogb.c
View File

@ -25,30 +25,39 @@ static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $";
#include <math.h>
#include <math_private.h>
int __ieee754_ilogb(double x)
int
__ieee754_ilogb (double x)
{
int32_t hx,lx,ix;
int32_t hx, lx, ix;
GET_HIGH_WORD(hx,x);
hx &= 0x7fffffff;
if(hx<0x00100000) {
GET_LOW_WORD(lx,x);
if((hx|lx)==0)
return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
else /* subnormal x */
if(hx==0) {
for (ix = -1043; lx>0; lx<<=1) ix -=1;
} else {
for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
}
return ix;
GET_HIGH_WORD (hx, x);
hx &= 0x7fffffff;
if (hx < 0x00100000)
{
GET_LOW_WORD (lx, x);
if ((hx | lx) == 0)
return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
else /* subnormal x */
if (hx == 0)
{
for (ix = -1043; lx > 0; lx <<= 1)
ix -= 1;
}
else if (hx<0x7ff00000) return (hx>>20)-1023;
else if (FP_ILOGBNAN != INT_MAX) {
/* ISO C99 requires ilogb(+-Inf) == INT_MAX. */
GET_LOW_WORD(lx,x);
if(((hx^0x7ff00000)|lx) == 0)
return INT_MAX;
else
{
for (ix = -1022, hx <<= 11; hx > 0; hx <<= 1)
ix -= 1;
}
return FP_ILOGBNAN;
return ix;
}
else if (hx < 0x7ff00000)
return (hx >> 20) - 1023;
else if (FP_ILOGBNAN != INT_MAX)
{
/* ISO C99 requires ilogb(+-Inf) == INT_MAX. */
GET_LOW_WORD (lx, x);
if (((hx ^ 0x7ff00000) | lx) == 0)
return INT_MAX;
}
return FP_ILOGBNAN;
}

374
sysdeps/ieee754/dbl-64/e_j0.c
View File

@ -11,7 +11,7 @@
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
for performance improvement on pipelined processors.
*/
*/
/* __ieee754_j0(x), __ieee754_y0(x)
* Bessel function of the first and second kinds of order zero.
@ -61,144 +61,166 @@
#include <math.h>
#include <math_private.h>
static double pzero(double), qzero(double);
static double pzero (double), qzero (double);
static const double
huge = 1e300,
one = 1.0,
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0, 2.00] */
R[] = {0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
-1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
-4.61832688532103189199e-09}, /* 0xBE33D5E7, 0x73D63FCE */
S[] = {0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
1.16614003333790000205e-09}; /* 0x3E1408BC, 0xF4745D8F */
huge = 1e300,
one = 1.0,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0, 2.00] */
R[] = { 0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
-1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
-4.61832688532103189199e-09 }, /* 0xBE33D5E7, 0x73D63FCE */
S[] = { 0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
1.16614003333790000205e-09 }; /* 0x3E1408BC, 0xF4745D8F */
static const double zero = 0.0;
double
__ieee754_j0(double x)
__ieee754_j0 (double x)
{
double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4;
int32_t hx,ix;
double z, s, c, ss, cc, r, u, v, r1, r2, s1, s2, z2, z4;
int32_t hx, ix;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return one/(x*x);
x = fabs(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
__sincos (x, &s, &c);
ss = s-c;
cc = s+c;
if(ix<0x7fe00000) { /* make sure x+x not overflow */
z = -__cos(x+x);
if ((s*c)<zero) cc = z/ss;
else ss = z/cc;
}
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(x);
else {
u = pzero(x); v = qzero(x);
z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(x);
}
return z;
GET_HIGH_WORD (hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
return one / (x * x);
x = fabs (x);
if (ix >= 0x40000000) /* |x| >= 2.0 */
{
__sincos (x, &s, &c);
ss = s - c;
cc = s + c;
if (ix < 0x7fe00000) /* make sure x+x not overflow */
{
z = -__cos (x + x);
if ((s * c) < zero)
cc = z / ss;
else
ss = z / cc;
}
if(ix<0x3f200000) { /* |x| < 2**-13 */
math_force_eval(huge+x); /* raise inexact if x != 0 */
if(ix<0x3e400000) return one; /* |x|<2**-27 */
else return one - 0.25*x*x;
}
z = x*x;
r1 = z*R[2]; z2=z*z;
r2 = R[3]+z*R[4]; z4=z2*z2;
r = r1 + z2*r2 + z4*R[5];
s1 = one+z*S[1];
s2 = S[2]+z*S[3];
s = s1 + z2*s2 + z4*S[4];
if(ix < 0x3FF00000) { /* |x| < 1.00 */
return one + z*(-0.25+(r/s));
} else {
u = 0.5*x;
return((one+u)*(one-u)+z*(r/s));
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix > 0x48000000)
z = (invsqrtpi * cc) / __ieee754_sqrt (x);
else
{
u = pzero (x); v = qzero (x);
z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (x);
}
return z;
}
if (ix < 0x3f200000) /* |x| < 2**-13 */
{
math_force_eval (huge + x); /* raise inexact if x != 0 */
if (ix < 0x3e400000)
return one; /* |x|<2**-27 */
else
return one - 0.25 * x * x;
}
z = x * x;
r1 = z * R[2]; z2 = z * z;
r2 = R[3] + z * R[4]; z4 = z2 * z2;
r = r1 + z2 * r2 + z4 * R[5];
s1 = one + z * S[1];
s2 = S[2] + z * S[3];
s = s1 + z2 * s2 + z4 * S[4];
if (ix < 0x3FF00000) /* |x| < 1.00 */
{
return one + z * (-0.25 + (r / s));
}
else
{
u = 0.5 * x;
return ((one + u) * (one - u) + z * (r / s));
}
}
strong_alias (__ieee754_j0, __j0_finite)
static const double
U[] = {-7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
-1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
-3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
-3.98205194132103398453e-11}, /* 0xBDC5E43D, 0x693FB3C8 */
V[] = {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
4.41110311332675467403e-10}; /* 0x3DFE5018, 0x3BD6D9EF */
U[] = { -7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
-1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
-3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
-3.98205194132103398453e-11 }, /* 0xBDC5E43D, 0x693FB3C8 */
V[] = { 1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
4.41110311332675467403e-10 }; /* 0x3DFE5018, 0x3BD6D9EF */
double
__ieee754_y0(double x)
__ieee754_y0 (double x)
{
double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2;
int32_t hx,ix,lx;
double z, s, c, ss, cc, u, v, z2, z4, z6, u1, u2, u3, v1, v2;
int32_t hx, ix, lx;
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */
if(hx<0) return zero/(zero*x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) + cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
__sincos (x, &s, &c);
ss = s-c;
cc = s+c;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if(ix<0x7fe00000) { /* make sure x+x not overflow */
z = -__cos(x+x);
if ((s*c)<zero) cc = z/ss;
else ss = z/cc;
}
if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
else {
u = pzero(x); v = qzero(x);
z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
}
return z;
EXTRACT_WORDS (hx, lx, x);
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
if (ix >= 0x7ff00000)
return one / (x + x * x);
if ((ix | lx) == 0)
return -HUGE_VAL + x; /* -inf and overflow exception. */
if (hx < 0)
return zero / (zero * x);
if (ix >= 0x40000000) /* |x| >= 2.0 */
{ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) + cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
__sincos (x, &s, &c);
ss = s - c;
cc = s + c;
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix < 0x7fe00000) /* make sure x+x not overflow */
{
z = -__cos (x + x);
if ((s * c) < zero)
cc = z / ss;
else
ss = z / cc;
}
if(ix<=0x3e400000) { /* x < 2**-27 */
return(U[0] + tpi*__ieee754_log(x));
if (ix > 0x48000000)
z = (invsqrtpi * ss) / __ieee754_sqrt (x);
else
{
u = pzero (x); v = qzero (x);
z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x);
}
z = x*x;
u1 = U[0]+z*U[1]; z2=z*z;
u2 = U[2]+z*U[3]; z4=z2*z2;
u3 = U[4]+z*U[5]; z6=z4*z2;
u = u1 + z2*u2 + z4*u3 + z6*U[6];
v1 = one+z*V[0];
v2 = V[1]+z*V[2];
v = v1 + z2*v2 + z4*V[3];
return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
return z;
}
if (ix <= 0x3e400000) /* x < 2**-27 */
{
return (U[0] + tpi * __ieee754_log (x));
}
z = x * x;
u1 = U[0] + z * U[1]; z2 = z * z;
u2 = U[2] + z * U[3]; z4 = z2 * z2;
u3 = U[4] + z * U[5]; z6 = z4 * z2;
u = u1 + z2 * u2 + z4 * u3 + z6 * U[6];
v1 = one + z * V[0];
v2 = V[1] + z * V[2];
v = v1 + z2 * v2 + z4 * V[3];
return (u / v + tpi * (__ieee754_j0 (x) * __ieee754_log (x)));
}
strong_alias (__ieee754_y0, __y0_finite)
@ -276,28 +298,43 @@ static const double pS2[5] = {
};
static double
pzero(double x)
pzero (double x)
{
const double *p,*q;
double z,r,s,z2,z4,r1,r2,r3,s1,s2,s3;
int32_t ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
if (ix>=0x41b00000) {return one;}
else if(ix>=0x40200000){p = pR8; q= pS8;}
else if(ix>=0x40122E8B){p = pR5; q= pS5;}
else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
else if(ix>=0x40000000){p = pR2; q= pS2;}
z = one/(x*x);
r1 = p[0]+z*p[1]; z2=z*z;
r2 = p[2]+z*p[3]; z4=z2*z2;
r3 = p[4]+z*p[5];
r = r1 + z2*r2 + z4*r3;
s1 = one+z*q[0];
s2 = q[1]+z*q[2];
s3 = q[3]+z*q[4];
s = s1 + z2*s2 + z4*s3;
return one+ r/s;
const double *p, *q;
double z, r, s, z2, z4, r1, r2, r3, s1, s2, s3;
int32_t ix;
GET_HIGH_WORD (ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41b00000)
{
return one;
}
else if (ix >= 0x40200000)
{
p = pR8; q = pS8;
}
else if (ix >= 0x40122E8B)
{
p = pR5; q = pS5;
}
else if (ix >= 0x4006DB6D)
{
p = pR3; q = pS3;
}
else if (ix >= 0x40000000)
{
p = pR2; q = pS2;
}
z = one / (x * x);
r1 = p[0] + z * p[1]; z2 = z * z;
r2 = p[2] + z * p[3]; z4 = z2 * z2;
r3 = p[4] + z * p[5];
r = r1 + z2 * r2 + z4 * r3;
s1 = one + z * q[0];
s2 = q[1] + z * q[2];
s3 = q[3] + z * q[4];
s = s1 + z2 * s2 + z4 * s3;
return one + r / s;
}
@ -379,26 +416,41 @@ static const double qS2[6] = {
};
static double
qzero(double x)
qzero (double x)
{
const double *p,*q;
double s,r,z,z2,z4,z6,r1,r2,r3,s1,s2,s3;
int32_t ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
if (ix>=0x41b00000) {return -.125/x;}
else if(ix>=0x40200000){p = qR8; q= qS8;}
else if(ix>=0x40122E8B){p = qR5; q= qS5;}
else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
else if(ix>=0x40000000){p = qR2; q= qS2;}
z = one/(x*x);
r1 = p[0]+z*p[1]; z2=z*z;
r2 = p[2]+z*p[3]; z4=z2*z2;
r3 = p[4]+z*p[5]; z6=z4*z2;
r= r1 + z2*r2 + z4*r3;
s1 = one+z*q[0];
s2 = q[1]+z*q[2];
s3 = q[3]+z*q[4];
s = s1 + z2*s2 + z4*s3 +z6*q[5];
return (-.125 + r/s)/x;
const double *p, *q;
double s, r, z, z2, z4, z6, r1, r2, r3, s1, s2, s3;
int32_t ix;
GET_HIGH_WORD (ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41b00000)
{
return -.125 / x;
}
else if (ix >= 0x40200000)
{
p = qR8; q = qS8;
}
else if (ix >= 0x40122E8B)
{
p = qR5; q = qS5;
}
else if (ix >= 0x4006DB6D)
{
p = qR3; q = qS3;
}
else if (ix >= 0x40000000)
{
p = qR2; q = qS2;
}
z = one / (x * x);
r1 = p[0] + z * p[1]; z2 = z * z;
r2 = p[2] + z * p[3]; z4 = z2 * z2;
r3 = p[4] + z * p[5]; z6 = z4 * z2;
r = r1 + z2 * r2 + z4 * r3;
s1 = one + z * q[0];
s2 = q[1] + z * q[2];
s3 = q[3] + z * q[4];
s = s1 + z2 * s2 + z4 * s3 + z6 * q[5];
return (-.125 + r / s) / x;
}

345
sysdeps/ieee754/dbl-64/e_j1.c
View File

@ -11,7 +11,7 @@
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
for performance improvement on pipelined processors.
*/
*/
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
@ -61,70 +61,81 @@
#include <math.h>
#include <math_private.h>
static double pone(double), qone(double);
static double pone (double), qone (double);
static const double
huge = 1e300,
one = 1.0,
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0,2] */
R[] = {-6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
-1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
4.96727999609584448412e-08}, /* 0x3E6AAAFA, 0x46CA0BD9 */
S[] = {0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
1.23542274426137913908e-11}; /* 0x3DAB2ACF, 0xCFB97ED8 */
huge = 1e300,
one = 1.0,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0,2] */
R[] = { -6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */
1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */
-1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */
4.96727999609584448412e-08 }, /* 0x3E6AAAFA, 0x46CA0BD9 */
S[] = { 0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */
1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
1.23542274426137913908e-11 }; /* 0x3DAB2ACF, 0xCFB97ED8 */
static const double zero = 0.0;
static const double zero = 0.0;
double
__ieee754_j1(double x)
__ieee754_j1 (double x)
{
double z, s,c,ss,cc,r,u,v,y,r1,r2,s1,s2,s3,z2,z4;
int32_t hx,ix;
double z, s, c, ss, cc, r, u, v, y, r1, r2, s1, s2, s3, z2, z4;
int32_t hx, ix;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
if(__builtin_expect(ix>=0x7ff00000, 0)) return one/x;
y = fabs(x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
__sincos (y, &s, &c);
ss = -s-c;
cc = s-c;
if(ix<0x7fe00000) { /* make sure y+y not overflow */
z = __cos(y+y);
if ((s*c)>zero) cc = z/ss;
else ss = z/cc;
}
/*
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
*/
if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(y);
else {
u = pone(y); v = qone(y);
z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(y);
}
if(hx<0) return -z;
else return z;
GET_HIGH_WORD (hx, x);
ix = hx & 0x7fffffff;
if (__builtin_expect (ix >= 0x7ff00000, 0))
return one / x;
y = fabs (x);
if (ix >= 0x40000000) /* |x| >= 2.0 */
{
__sincos (y, &s, &c);
ss = -s - c;
cc = s - c;
if (ix < 0x7fe00000) /* make sure y+y not overflow */
{
z = __cos (y + y);
if ((s * c) > zero)
cc = z / ss;
else
ss = z / cc;
}
if(__builtin_expect(ix<0x3e400000, 0)) { /* |x|<2**-27 */
if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
/*
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
*/
if (ix > 0x48000000)
z = (invsqrtpi * cc) / __ieee754_sqrt (y);
else
{
u = pone (y); v = qone (y);
z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (y);
}
z = x*x;
r1 = z*R[0]; z2=z*z;
r2 = R[1]+z*R[2]; z4=z2*z2;
r = r1 + z2*r2 + z4*R[3];
r *= x;
s1 = one+z*S[1];
s2 = S[2]+z*S[3];
s3 = S[4]+z*S[5];
s = s1 + z2*s2 + z4*s3;
return(x*0.5+r/s);
if (hx < 0)
return -z;
else
return z;
}
if (__builtin_expect (ix < 0x3e400000, 0)) /* |x|<2**-27 */
{
if (huge + x > one)
return 0.5 * x; /* inexact if x!=0 necessary */
}
z = x * x;
r1 = z * R[0]; z2 = z * z;
r2 = R[1] + z * R[2]; z4 = z2 * z2;
r = r1 + z2 * r2 + z4 * R[3];
r *= x;
s1 = one + z * S[1];
s2 = S[2] + z * S[3];
s3 = S[4] + z * S[5];
s = s1 + z2 * s2 + z4 * s3;
return (x * 0.5 + r / s);
}
strong_alias (__ieee754_j1, __j1_finite)
@ -144,57 +155,67 @@ static const double V0[5] = {
};
double
__ieee754_y1(double x)
__ieee754_y1 (double x)
{
double z, s,c,ss,cc,u,v,u1,u2,v1,v2,v3,z2,z4;
int32_t hx,ix,lx;
double z, s, c, ss, cc, u, v, u1, u2, v1, v2, v3, z2, z4;
int32_t hx, ix, lx;
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(__builtin_expect(ix>=0x7ff00000, 0)) return one/(x+x*x);
if(__builtin_expect((ix|lx)==0, 0))
return -HUGE_VAL+x; /* -inf and overflow exception. */;
if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */
__sincos (x, &s, &c);
ss = -s-c;
cc = s-c;
if(ix<0x7fe00000) { /* make sure x+x not overflow */
z = __cos(x+x);
if ((s*c)>zero) cc = z/ss;
else ss = z/cc;
}
/* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
* where x0 = x-3pi/4
* Better formula:
* cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (cos(x) + sin(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
else {
u = pone(x); v = qone(x);
z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
}
return z;
EXTRACT_WORDS (hx, lx, x);
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (__builtin_expect (ix >= 0x7ff00000, 0))
return one / (x + x * x);
if (__builtin_expect ((ix | lx) == 0, 0))
return -HUGE_VAL + x;
/* -inf and overflow exception. */;
if (__builtin_expect (hx < 0, 0))
return zero / (zero * x);
if (ix >= 0x40000000) /* |x| >= 2.0 */
{
__sincos (x, &s, &c);
ss = -s - c;
cc = s - c;
if (ix < 0x7fe00000) /* make sure x+x not overflow */
{
z = __cos (x + x);
if ((s * c) > zero)
cc = z / ss;
else
ss = z / cc;
}
if(__builtin_expect(ix<=0x3c900000, 0)) { /* x < 2**-54 */
return(-tpi/x);
/* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
* where x0 = x-3pi/4
* Better formula:
* cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
* = -1/sqrt(2) * (cos(x) + sin(x))
* To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
if (ix > 0x48000000)
z = (invsqrtpi * ss) / __ieee754_sqrt (x);
else
{
u = pone (x); v = qone (x);
z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x);
}
z = x*x;
u1 = U0[0]+z*U0[1];z2=z*z;
u2 = U0[2]+z*U0[3];z4=z2*z2;
u = u1 + z2*u2 + z4*U0[4];
v1 = one+z*V0[0];
v2 = V0[1]+z*V0[2];
v3 = V0[3]+z*V0[4];
v = v1 + z2*v2 + z4*v3;
return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
return z;
}
if (__builtin_expect (ix <= 0x3c900000, 0)) /* x < 2**-54 */
{
return (-tpi / x);
}
z = x * x;
u1 = U0[0] + z * U0[1]; z2 = z * z;
u2 = U0[2] + z * U0[3]; z4 = z2 * z2;
u = u1 + z2 * u2 + z4 * U0[4];
v1 = one + z * V0[0];
v2 = V0[1] + z * V0[2];
v3 = V0[3] + z * V0[4];
v = v1 + z2 * v2 + z4 * v3;
return (x * (u / v) + tpi * (__ieee754_j1 (x) * __ieee754_log (x) - one / x));
}
strong_alias (__ieee754_y1, __y1_finite)
@ -256,7 +277,7 @@ static const double ps3[5] = {
1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
};
static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
static const double pr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */
1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
@ -273,28 +294,43 @@ static const double ps2[5] = {
};
static double
pone(double x)
pone (double x)
{
const double *p,*q;
double z,r,s,r1,r2,r3,s1,s2,s3,z2,z4;
int32_t ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
if (ix>=0x41b00000) {return one;}
else if(ix>=0x40200000){p = pr8; q= ps8;}
else if(ix>=0x40122E8B){p = pr5; q= ps5;}
else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
else if(ix>=0x40000000){p = pr2; q= ps2;}
z = one/(x*x);
r1 = p[0]+z*p[1]; z2=z*z;
r2 = p[2]+z*p[3]; z4=z2*z2;
r3 = p[4]+z*p[5];
r = r1 + z2*r2 + z4*r3;
s1 = one+z*q[0];
s2 = q[1]+z*q[2];
s3 = q[3]+z*q[4];
s = s1 + z2*s2 + z4*s3;
return one+ r/s;
const double *p, *q;
double z, r, s, r1, r2, r3, s1, s2, s3, z2, z4;
int32_t ix;
GET_HIGH_WORD (ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41b00000)
{
return one;
}
else if (ix >= 0x40200000)
{
p = pr8; q = ps8;
}
else if (ix >= 0x40122E8B)
{
p = pr5; q = ps5;
}
else if (ix >= 0x4006DB6D)
{
p = pr3; q = ps3;
}
else if (ix >= 0x40000000)
{
p = pr2; q = ps2;
}
z = one / (x * x);
r1 = p[0] + z * p[1]; z2 = z * z;
r2 = p[2] + z * p[3]; z4 = z2 * z2;
r3 = p[4] + z * p[5];
r = r1 + z2 * r2 + z4 * r3;
s1 = one + z * q[0];
s2 = q[1] + z * q[2];
s3 = q[3] + z * q[4];
s = s1 + z2 * s2 + z4 * s3;
return one + r / s;
}
@ -359,7 +395,7 @@ static const double qs3[6] = {
-1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
};
static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
static const double qr2[6] = { /* for x in [2.8570,2]=1/[0.3499,0.5] */
-1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
-1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
-2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
@ -377,26 +413,41 @@ static const double qs2[6] = {
};
static double
qone(double x)
qone (double x)
{
const double *p,*q;
double s,r,z,r1,r2,r3,s1,s2,s3,z2,z4,z6;
int32_t ix;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
if (ix>=0x41b00000) {return .375/x;}
else if(ix>=0x40200000){p = qr8; q= qs8;}
else if(ix>=0x40122E8B){p = qr5; q= qs5;}
else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
else if(ix>=0x40000000){p = qr2; q= qs2;}
z = one/(x*x);
r1 = p[0]+z*p[1]; z2=z*z;
r2 = p[2]+z*p[3]; z4=z2*z2;
r3 = p[4]+z*p[5]; z6=z4*z2;
r = r1 + z2*r2 + z4*r3;
s1 = one+z*q[0];
s2 = q[1]+z*q[2];
s3 = q[3]+z*q[4];
s = s1 + z2*s2 + z4*s3 + z6*q[5];
return (.375 + r/s)/x;
const double *p, *q;
double s, r, z, r1, r2, r3, s1, s2, s3, z2, z4, z6;
int32_t ix;
GET_HIGH_WORD (ix, x);
ix &= 0x7fffffff;
if (ix >= 0x41b00000)
{
return .375 / x;
}
else if (ix >= 0x40200000)
{
p = qr8; q = qs8;
}
else if (ix >= 0x40122E8B)
{
p = qr5; q = qs5;
}
else if (ix >= 0x4006DB6D)
{
p = qr3; q = qs3;
}
else if (ix >= 0x40000000)
{
p = qr2; q = qs2;
}
z = one / (x * x);
r1 = p[0] + z * p[1]; z2 = z * z;
r2 = p[2] + z * p[3]; z4 = z2 * z2;
r3 = p[4] + z * p[5]; z6 = z4 * z2;
r = r1 + z2 * r2 + z4 * r3;
s1 = one + z * q[0];
s2 = q[1] + z * q[2];
s3 = q[3] + z * q[4];
s = s1 + z2 * s2 + z4 * s3 + z6 * q[5];
return (.375 + r / s) / x;
}

488
sysdeps/ieee754/dbl-64/e_jn.c
View File

@ -41,246 +41,284 @@
#include <math_private.h>
static const double
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
static const double zero = 0.00000000000000000000e+00;
static const double zero = 0.00000000000000000000e+00;
double
__ieee754_jn(int n, double x)
__ieee754_jn (int n, double x)
{
int32_t i,hx,ix,lx, sgn;
double a, b, temp, di;
double z, w;
int32_t i, hx, ix, lx, sgn;
double a, b, temp, di;
double z, w;
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if J(n,NaN) is NaN */
if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000, 0))
return x+x;
if(n<0){
n = -n;
x = -x;
hx ^= 0x80000000;
}
if(n==0) return(__ieee754_j0(x));
if(n==1) return(__ieee754_j1(x));
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabs(x);
if(__builtin_expect((ix|lx)==0||ix>=0x7ff00000,0))
/* if x is 0 or inf */
b = zero;
else if((double)n<=x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if(ix>=0x52D00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
double s;
double c;
__sincos (x, &s, &c);
switch(n&3) {
case 0: temp = c + s; break;
case 1: temp = -c + s; break;
case 2: temp = -c - s; break;
case 3: temp = c - s; break;
}
b = invsqrtpi*temp/__ieee754_sqrt(x);
} else {
a = __ieee754_j0(x);
b = __ieee754_j1(x);
for(i=1;i<n;i++){
temp = b;
b = b*((double)(i+i)/x) - a; /* avoid underflow */
a = temp;
}
/* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
* Thus, J(-n,x) = J(n,-x)
*/
EXTRACT_WORDS (hx, lx, x);
ix = 0x7fffffff & hx;
/* if J(n,NaN) is NaN */
if (__builtin_expect ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000, 0))
return x + x;
if (n < 0)
{
n = -n;
x = -x;
hx ^= 0x80000000;
}
if (n == 0)
return (__ieee754_j0 (x));
if (n == 1)
return (__ieee754_j1 (x));
sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
x = fabs (x);
if (__builtin_expect ((ix | lx) == 0 || ix >= 0x7ff00000, 0))
/* if x is 0 or inf */
b = zero;
else if ((double) n <= x)
{
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
if (ix >= 0x52D00000) /* x > 2**302 */
{ /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
double s;
double c;
__sincos (x, &s, &c);
switch (n & 3)
{
case 0: temp = c + s; break;
case 1: temp = -c + s; break;
case 2: temp = -c - s; break;
case 3: temp = c - s; break;
}
} else {
if(ix<0x3e100000) { /* x < 2**-29 */
/* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if(n>33) /* underflow */
b = zero;
else {
temp = x*0.5; b = temp;
for (a=one,i=2;i<=n;i++) {
a *= (double)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
b = b/a;
}
} else {
/* use backward recurrence */
/* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
* 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
* -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction:
* 1
* = -----------------------
* 1
* w - -----------------
* 1
* w+h - ---------
* w+2h - ...
*
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple
*/
/* determine k */
double t,v;
double q0,q1,h,tmp; int32_t k,m;
w = (n+n)/(double)x; h = 2.0/(double)x;
q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
while(q1<1.0e9) {
k += 1; z += h;
tmp = z*q1 - q0;
q0 = q1;
q1 = tmp;
}
m = n+n;
for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
a = t;
b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = n;
v = two/x;
tmp = tmp*__ieee754_log(fabs(v*tmp));
if(tmp<7.09782712893383973096e+02) {
for(i=n-1,di=(double)(i+i);i>0;i--){
temp = b;
b *= di;
b = b/x - a;
a = temp;
di -= two;
}
} else {
for(i=n-1,di=(double)(i+i);i>0;i--){
temp = b;
b *= di;
b = b/x - a;
a = temp;
di -= two;
/* scale b to avoid spurious overflow */
if(b>1e100) {
a /= b;
t /= b;
b = one;
}
}
}
/* j0() and j1() suffer enormous loss of precision at and
* near zero; however, we know that their zero points never
* coincide, so just choose the one further away from zero.
*/
z = __ieee754_j0 (x);
w = __ieee754_j1 (x);
if (fabs (z) >= fabs (w))
b = (t * z / b);
else
b = (t * w / a);
b = invsqrtpi * temp / __ieee754_sqrt (x);
}
else
{
a = __ieee754_j0 (x);
b = __ieee754_j1 (x);
for (i = 1; i < n; i++)
{
temp = b;
b = b * ((double) (i + i) / x) - a; /* avoid underflow */
a = temp;
}
}
if(sgn==1) return -b; else return b;
}
else
{
if (ix < 0x3e100000) /* x < 2**-29 */
{ /* x is tiny, return the first Taylor expansion of J(n,x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
b = zero;
else
{
temp = x * 0.5; b = temp;
for (a = one, i = 2; i <= n; i++)
{
a *= (double) i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
b = b / a;
}
}
else
{
/* use backward recurrence */
/* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
* 2n - 2(n+1) - 2(n+2)
*
* 1 1 1
* (for large x) = ---- ------ ------ .....
* 2n 2(n+1) 2(n+2)
* -- - ------ - ------ -
* x x x
*
* Let w = 2n/x and h=2/x, then the above quotient
* is equal to the continued fraction:
* 1
* = -----------------------
* 1
* w - -----------------
* 1
* w+h - ---------
* w+2h - ...
*
* To determine how many terms needed, let
* Q(0) = w, Q(1) = w(w+h) - 1,
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
* When Q(k) > 1e4 good for single
* When Q(k) > 1e9 good for double
* When Q(k) > 1e17 good for quadruple
*/
/* determine k */
double t, v;
double q0, q1, h, tmp; int32_t k, m;
w = (n + n) / (double) x; h = 2.0 / (double) x;
q0 = w; z = w + h; q1 = w * z - 1.0; k = 1;
while (q1 < 1.0e9)
{
k += 1; z += h;
tmp = z * q1 - q0;
q0 = q1;
q1 = tmp;
}
m = n + n;
for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
t = one / (i / x - t);
a = t;
b = one;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* long double 1.1356523406294143949491931077970765006170e+04
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
tmp = n;
v = two / x;
tmp = tmp * __ieee754_log (fabs (v * tmp));
if (tmp < 7.09782712893383973096e+02)
{
for (i = n - 1, di = (double) (i + i); i > 0; i--)
{
temp = b;
b *= di;
b = b / x - a;
a = temp;
di -= two;
}
}
else
{
for (i = n - 1, di = (double) (i + i); i > 0; i--)
{
temp = b;
b *= di;
b = b / x - a;
a = temp;
di -= two;
/* scale b to avoid spurious overflow */
if (b > 1e100)
{
a /= b;
t /= b;
b = one;
}
}
}
/* j0() and j1() suffer enormous loss of precision at and
* near zero; however, we know that their zero points never
* coincide, so just choose the one further away from zero.
*/
z = __ieee754_j0 (x);
w = __ieee754_j1 (x);
if (fabs (z) >= fabs (w))
b = (t * z / b);
else
b = (t * w / a);
}
}
if (sgn == 1)
return -b;
else
return b;
}
strong_alias (__ieee754_jn, __jn_finite)
double
__ieee754_yn(int n, double x)
__ieee754_yn (int n, double x)
{
int32_t i,hx,ix,lx;
int32_t sign;
double a, b, temp;
int32_t i, hx, ix, lx;
int32_t sign;
double a, b, temp;
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */
if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000,0))
return x+x;
if(__builtin_expect((ix|lx)==0, 0))
return -HUGE_VAL+x; /* -inf and overflow exception. */;
if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
sign = 1;
if(n<0){
n = -n;
sign = 1 - ((n&1)<<1);
EXTRACT_WORDS (hx, lx, x);
ix = 0x7fffffff & hx;
/* if Y(n,NaN) is NaN */
if (__builtin_expect ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000, 0))
return x + x;
if (__builtin_expect ((ix | lx) == 0, 0))
return -HUGE_VAL + x;
/* -inf and overflow exception. */;
if (__builtin_expect (hx < 0, 0))
return zero / (zero * x);
sign = 1;
if (n < 0)
{
n = -n;
sign = 1 - ((n & 1) << 1);
}
if (n == 0)
return (__ieee754_y0 (x));
if (n == 1)
return (sign * __ieee754_y1 (x));
if (__builtin_expect (ix == 0x7ff00000, 0))
return zero;
if (ix >= 0x52D00000) /* x > 2**302 */
{ /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
double c;
double s;
__sincos (x, &s, &c);
switch (n & 3)
{
case 0: temp = s - c; break;
case 1: temp = -s - c; break;
case 2: temp = -s + c; break;
case 3: temp = s + c; break;
}
if(n==0) return(__ieee754_y0(x));
if(n==1) return(sign*__ieee754_y1(x));
if(__builtin_expect(ix==0x7ff00000, 0)) return zero;
if(ix>=0x52D00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
*
* n sin(xn)*sqt2 cos(xn)*sqt2
* ----------------------------------
* 0 s-c c+s
* 1 -s-c -c+s
* 2 -s+c -c-s
* 3 s+c c-s
*/
double c;
double s;
__sincos (x, &s, &c);
switch(n&3) {
case 0: temp = s - c; break;
case 1: temp = -s - c; break;
case 2: temp = -s + c; break;
case 3: temp = s + c; break;
}
b = invsqrtpi*temp/__ieee754_sqrt(x);
} else {
u_int32_t high;
a = __ieee754_y0(x);
b = __ieee754_y1(x);
/* quit if b is -inf */
GET_HIGH_WORD(high,b);
for(i=1;i<n&&high!=0xfff00000;i++){
temp = b;
b = ((double)(i+i)/x)*b - a;
GET_HIGH_WORD(high,b);
a = temp;
}
/* If B is +-Inf, set up errno accordingly. */
if (! __finite (b))
__set_errno (ERANGE);
b = invsqrtpi * temp / __ieee754_sqrt (x);
}
else
{
u_int32_t high;
a = __ieee754_y0 (x);
b = __ieee754_y1 (x);
/* quit if b is -inf */
GET_HIGH_WORD (high, b);
for (i = 1; i < n && high != 0xfff00000; i++)
{
temp = b;
b = ((double) (i + i) / x) * b - a;
GET_HIGH_WORD (high, b);
a = temp;
}
if(sign>0) return b; else return -b;
/* If B is +-Inf, set up errno accordingly. */
if (!__finite (b))
__set_errno (ERANGE);
}
if (sign > 0)
return b;
else
return -b;
}
strong_alias (__ieee754_yn, __yn_finite)

16
sysdeps/ieee754/dbl-64/e_log.c
View File

@ -56,12 +56,12 @@ SECTION
__ieee754_log (double x)
{
#define M 4
static const int pr[M] = {8, 10, 18, 32};
static const int pr[M] = { 8, 10, 18, 32 };
int i, j, n, ux, dx, p;
double dbl_n, u, p0, q, r0, w, nln2a, luai, lubi, lvaj, lvbj,
sij, ssij, ttij, A, B, B0, y, y1, y2, polI, polII, sa, sb,
t1, t2, t7, t8, t, ra, rb, ww,
a0, aa0, s1, s2, ss2, s3, ss3, a1, aa1, a, aa, b, bb, c;
sij, ssij, ttij, A, B, B0, y, y1, y2, polI, polII, sa, sb,
t1, t2, t7, t8, t, ra, rb, ww,
a0, aa0, s1, s2, ss2, s3, ss3, a1, aa1, a, aa, b, bb, c;
#ifndef DLA_FMS
double t3, t4, t5, t6;
#endif
@ -80,15 +80,15 @@ __ieee754_log (double x)
if (__builtin_expect (ux < 0x00100000, 0))
{
if (__builtin_expect (((ux & 0x7fffffff) | dx) == 0, 0))
return MHALF / 0.0; /* return -INF */
return MHALF / 0.0; /* return -INF */
if (__builtin_expect (ux < 0, 0))
return (x - x) / 0.0; /* return NaN */
return (x - x) / 0.0; /* return NaN */
n -= 54;
x *= two54.d; /* scale x */
x *= two54.d; /* scale x */
num.d = x;
}
if (__builtin_expect (ux >= 0x7ff00000, 0))
return x + x; /* INF or NaN */
return x + x; /* INF or NaN */
/* Regular values of x */

16
sysdeps/ieee754/dbl-64/e_log10.c
View File

@ -46,10 +46,10 @@
#include <math.h>
#include <math_private.h>
static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
static const double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B, 0x1526E50E */
static const double log10_2hi = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
static const double log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
static const double ivln10 = 4.34294481903251816668e-01; /* 0x3FDBCB7B, 0x1526E50E */
static const double log10_2hi = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
static const double log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
double
__ieee754_log10 (double x)
@ -62,13 +62,13 @@ __ieee754_log10 (double x)
k = 0;
if (hx < 0x00100000)
{ /* x < 2**-1022 */
{ /* x < 2**-1022 */
if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0))
return -two54 / (x - x); /* log(+-0)=-inf */
return -two54 / (x - x); /* log(+-0)=-inf */
if (__builtin_expect (hx < 0, 0))
return (x - x) / (x - x); /* log(-#) = NaN */
return (x - x) / (x - x); /* log(-#) = NaN */
k -= 54;
x *= two54; /* subnormal number, scale up x */
x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD (hx, x);
}
if (__builtin_expect (hx >= 0x7ff00000, 0))

28
sysdeps/ieee754/dbl-64/e_log2.c
View File

@ -58,14 +58,14 @@
#include <math_private.h>
static const double ln2 = 0.69314718055994530942;
static const double two54 = 1.80143985094819840000e+16; /* 43500000 00000000 */
static const double Lg1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
static const double Lg2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
static const double Lg3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
static const double Lg4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
static const double Lg5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
static const double Lg6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
static const double Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double two54 = 1.80143985094819840000e+16; /* 43500000 00000000 */
static const double Lg1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
static const double Lg2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
static const double Lg3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
static const double Lg4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
static const double Lg5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
static const double Lg6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
static const double Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
@ -80,13 +80,13 @@ __ieee754_log2 (double x)
k = 0;
if (hx < 0x00100000)
{ /* x < 2**-1022 */
{ /* x < 2**-1022 */
if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0))
return -two54 / (x - x); /* log(+-0)=-inf */
return -two54 / (x - x); /* log(+-0)=-inf */
if (__builtin_expect (hx < 0, 0))
return (x - x) / (x - x); /* log(-#) = NaN */
return (x - x) / (x - x); /* log(-#) = NaN */
k -= 54;
x *= two54; /* subnormal number, scale up x */
x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD (hx, x);
}
if (__builtin_expect (hx >= 0x7ff00000, 0))
@ -94,12 +94,12 @@ __ieee754_log2 (double x)
k += (hx >> 20) - 1023;
hx &= 0x000fffff;
i = (hx + 0x95f64) & 0x100000;
SET_HIGH_WORD (x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
SET_HIGH_WORD (x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
k += (i >> 20);
dk = (double) k;
f = x - 1.0;
if ((0x000fffff & (2 + hx)) < 3)
{ /* |f| < 2**-20 */
{ /* |f| < 2**-20 */
if (f == zero)
return dk;
R = f * f * (0.5 - 0.33333333333333333 * f);

219
sysdeps/ieee754/dbl-64/e_rem_pio2.c
View File

@ -57,110 +57,137 @@ static const int32_t npio2_hw[] = {
*/
static const double
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
int32_t
__ieee754_rem_pio2(double x, double *y)
__ieee754_rem_pio2 (double x, double *y)
{
double z,w,t,r,fn;
double tx[3];
int32_t e0,i,j,nx,n,ix,hx;
u_int32_t low;
double z, w, t, r, fn;
double tx[3];
int32_t e0, i, j, nx, n, ix, hx;
u_int32_t low;
GET_HIGH_WORD(hx,x); /* high word of x */
ix = hx&0x7fffffff;
if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{y[0] = x; y[1] = 0; return 0;}
if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
if(hx>0) {
z = x - pio2_1;
if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z - pio2_1t;
y[1] = (z-y[0])-pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z-y[0])-pio2_2t;
GET_HIGH_WORD (hx, x); /* high word of x */
ix = hx & 0x7fffffff;
if (ix <= 0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{
y[0] = x; y[1] = 0; return 0;
}
if (ix < 0x4002d97c) /* |x| < 3pi/4, special case with n=+-1 */
{
if (hx > 0)
{
z = x - pio2_1;
if (ix != 0x3ff921fb) /* 33+53 bit pi is good enough */
{
y[0] = z - pio2_1t;
y[1] = (z - y[0]) - pio2_1t;
}
else /* near pi/2, use 33+33+53 bit pi */
{
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z - y[0]) - pio2_2t;
}
return 1;
}
else /* negative x */
{
z = x + pio2_1;
if (ix != 0x3ff921fb) /* 33+53 bit pi is good enough */
{
y[0] = z + pio2_1t;
y[1] = (z - y[0]) + pio2_1t;
}
else /* near pi/2, use 33+33+53 bit pi */
{
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z - y[0]) + pio2_2t;
}
return -1;
}
}
if (ix <= 0x413921fb) /* |x| ~<= 2^19*(pi/2), medium size */
{
t = fabs (x);
n = (int32_t) (t * invpio2 + half);
fn = (double) n;
r = t - fn * pio2_1;
w = fn * pio2_1t; /* 1st round good to 85 bit */
if (n < 32 && ix != npio2_hw[n - 1])
{
y[0] = r - w; /* quick check no cancellation */
}
else
{
u_int32_t high;
j = ix >> 20;
y[0] = r - w;
GET_HIGH_WORD (high, y[0]);
i = j - ((high >> 20) & 0x7ff);
if (i > 16) /* 2nd iteration needed, good to 118 */
{
t = r;
w = fn * pio2_2;
r = t - w;
w = fn * pio2_2t - ((t - r) - w);
y[0] = r - w;
GET_HIGH_WORD (high, y[0]);
i = j - ((high >> 20) & 0x7ff);
if (i > 49) /* 3rd iteration need, 151 bits acc */
{
t = r; /* will cover all possible cases */
w = fn * pio2_3;
r = t - w;
w = fn * pio2_3t - ((t - r) - w);
y[0] = r - w;
}
return 1;
} else { /* negative x */
z = x + pio2_1;
if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z + pio2_1t;
y[1] = (z-y[0])+pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z-y[0])+pio2_2t;
}
return -1;
}
}
if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
t = fabs(x);
n = (int32_t) (t*invpio2+half);
fn = (double)n;
r = t-fn*pio2_1;
w = fn*pio2_1t; /* 1st round good to 85 bit */
if(n<32&&ix!=npio2_hw[n-1]) {
y[0] = r-w; /* quick check no cancellation */
} else {
u_int32_t high;
j = ix>>20;
y[0] = r-w;
GET_HIGH_WORD(high,y[0]);
i = j-((high>>20)&0x7ff);
if(i>16) { /* 2nd iteration needed, good to 118 */
t = r;
w = fn*pio2_2;
r = t-w;
w = fn*pio2_2t-((t-r)-w);
y[0] = r-w;
GET_HIGH_WORD(high,y[0]);
i = j-((high>>20)&0x7ff);
if(i>49) { /* 3rd iteration need, 151 bits acc */
t = r; /* will cover all possible cases */
w = fn*pio2_3;
r = t-w;
w = fn*pio2_3t-((t-r)-w);
y[0] = r-w;
}
}
}
y[1] = (r-y[0])-w;
if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
else return n;
y[1] = (r - y[0]) - w;
if (hx < 0)
{
y[0] = -y[0]; y[1] = -y[1]; return -n;
}
/*
* all other (large) arguments
*/
if(ix>=0x7ff00000) { /* x is inf or NaN */
y[0]=y[1]=x-x; return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
GET_LOW_WORD(low,x);
SET_LOW_WORD(z,low);
e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
for(i=0;i<2;i++) {
tx[i] = (double)((int32_t)(z));
z = (z-tx[i])*two24;
}
tx[2] = z;
nx = 3;
while(tx[nx-1]==zero) nx--; /* skip zero term */
n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
else
return n;
}
/*
* all other (large) arguments
*/
if (ix >= 0x7ff00000) /* x is inf or NaN */
{
y[0] = y[1] = x - x; return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
GET_LOW_WORD (low, x);
SET_LOW_WORD (z, low);
e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */
SET_HIGH_WORD (z, ix - ((int32_t) (e0 << 20)));
for (i = 0; i < 2; i++)
{
tx[i] = (double) ((int32_t) (z));
z = (z - tx[i]) * two24;
}
tx[2] = z;
nx = 3;
while (tx[nx - 1] == zero)
nx--; /* skip zero term */
n = __kernel_rem_pio2 (tx, y, e0, nx, 2, two_over_pi);
if (hx < 0)
{
y[0] = -y[0]; y[1] = -y[1]; return -n;
}
return n;
}
#endif

182
sysdeps/ieee754/dbl-64/e_remainder.c
View File

@ -39,89 +39,111 @@
/* An ultimate remainder routine. Given two IEEE double machine numbers x */
/* ,y it computes the correctly rounded (to nearest) value of remainder */
/**************************************************************************/
double __ieee754_remainder(double x, double y)
double
__ieee754_remainder (double x, double y)
{
double z,d,xx;
int4 kx,ky,n,nn,n1,m1,l;
mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
u.x=x;
t.x=y;
kx=u.i[HIGH_HALF]&0x7fffffff; /* no sign for x*/
t.i[HIGH_HALF]&=0x7fffffff; /*no sign for y */
ky=t.i[HIGH_HALF];
double z, d, xx;
int4 kx, ky, n, nn, n1, m1, l;
mynumber u, t, w = { { 0, 0 } }, v = { { 0, 0 } }, ww = { { 0, 0 } }, r;
u.x = x;
t.x = y;
kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign for x*/
t.i[HIGH_HALF] &= 0x7fffffff; /*no sign for y */
ky = t.i[HIGH_HALF];
/*------ |x| < 2^1023 and 2^-970 < |y| < 2^1024 ------------------*/
if (kx<0x7fe00000 && ky<0x7ff00000 && ky>=0x03500000) {
SET_RESTORE_ROUND_NOEX (FE_TONEAREST);
if (kx+0x00100000<ky) return x;
if ((kx-0x01500000)<ky) {
z=x/t.x;
v.i[HIGH_HALF]=t.i[HIGH_HALF];
d=(z+big.x)-big.x;
xx=(x-d*v.x)-d*(t.x-v.x);
if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
else {
if (ABS(xx)>0.5*t.x) return (z>d)?xx-t.x:xx+t.x;
else return xx;
}
} /* (kx<(ky+0x01500000)) */
else {
r.x=1.0/t.x;
n=t.i[HIGH_HALF];
nn=(n&0x7ff00000)+0x01400000;
w.i[HIGH_HALF]=n;
ww.x=t.x-w.x;
l=(kx-nn)&0xfff00000;
n1=ww.i[HIGH_HALF];
m1=r.i[HIGH_HALF];
while (l>0) {
r.i[HIGH_HALF]=m1-l;
z=u.x*r.x;
w.i[HIGH_HALF]=n+l;
ww.i[HIGH_HALF]=(n1)?n1+l:n1;
d=(z+big.x)-big.x;
u.x=(u.x-d*w.x)-d*ww.x;
l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
}
r.i[HIGH_HALF]=m1;
w.i[HIGH_HALF]=n;
ww.i[HIGH_HALF]=n1;
z=u.x*r.x;
d=(z+big.x)-big.x;
u.x=(u.x-d*w.x)-d*ww.x;
if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
else
if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
else
{z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
}
} /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */
else {
if (kx<0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
y=ABS(y)*t128.x;
z=__ieee754_remainder(x,y)*t128.x;
z=__ieee754_remainder(z,y)*tm128.x;
return z;
}
else {
if ((kx&0x7ff00000)==0x7fe00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
y=ABS(y);
z=2.0*__ieee754_remainder(0.5*x,y);
d = ABS(z);
if (d <= ABS(d-y)) return z;
else return (z>0)?z-y:z+y;
}
else { /* if x is too big */
if (ky==0 && t.i[LOW_HALF] == 0) /* y = 0 */
return (x * y) / (x * y);
else if (kx >= 0x7ff00000 /* x not finite */
|| (ky>0x7ff00000 /* y is NaN */
|| (ky == 0x7ff00000 && t.i[LOW_HALF] != 0)))
return (x * y) / (x * y);
else
if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000)
{
SET_RESTORE_ROUND_NOEX (FE_TONEAREST);
if (kx + 0x00100000 < ky)
return x;
if ((kx - 0x01500000) < ky)
{
z = x / t.x;
v.i[HIGH_HALF] = t.i[HIGH_HALF];
d = (z + big.x) - big.x;
xx = (x - d * v.x) - d * (t.x - v.x);
if (d - z != 0.5 && d - z != -0.5)
return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x);
else
{
if (ABS (xx) > 0.5 * t.x)
return (z > d) ? xx - t.x : xx + t.x;
else
return xx;
}
} /* (kx<(ky+0x01500000)) */
else
{
r.x = 1.0 / t.x;
n = t.i[HIGH_HALF];
nn = (n & 0x7ff00000) + 0x01400000;
w.i[HIGH_HALF] = n;
ww.x = t.x - w.x;
l = (kx - nn) & 0xfff00000;
n1 = ww.i[HIGH_HALF];
m1 = r.i[HIGH_HALF];
while (l > 0)
{
r.i[HIGH_HALF] = m1 - l;
z = u.x * r.x;
w.i[HIGH_HALF] = n + l;
ww.i[HIGH_HALF] = (n1) ? n1 + l : n1;
d = (z + big.x) - big.x;
u.x = (u.x - d * w.x) - d * ww.x;
l = (u.i[HIGH_HALF] & 0x7ff00000) - nn;
}
r.i[HIGH_HALF] = m1;
w.i[HIGH_HALF] = n;
ww.i[HIGH_HALF] = n1;
z = u.x * r.x;
d = (z + big.x) - big.x;
u.x = (u.x - d * w.x) - d * ww.x;
if (ABS (u.x) < 0.5 * t.x)
return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x);
else
if (ABS (u.x) > 0.5 * t.x)
return (d > z) ? u.x + t.x : u.x - t.x;
else
{
z = u.x / t.x; d = (z + big.x) - big.x;
return ((u.x - d * w.x) - d * ww.x);
}
}
} /* (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000) */
else
{
if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0))
{
y = ABS (y) * t128.x;
z = __ieee754_remainder (x, y) * t128.x;
z = __ieee754_remainder (z, y) * tm128.x;
return z;
}
else
{
if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 &&
(ky > 0 || t.i[LOW_HALF] != 0))
{
y = ABS (y);
z = 2.0 * __ieee754_remainder (0.5 * x, y);
d = ABS (z);
if (d <= ABS (d - y))
return z;
else
return (z > 0) ? z - y : z + y;
}
else /* if x is too big */
{
if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */
return (x * y) / (x * y);
else if (kx >= 0x7ff00000 /* x not finite */
|| (ky > 0x7ff00000 /* y is NaN */
|| (ky == 0x7ff00000 && t.i[LOW_HALF] != 0)))
return (x * y) / (x * y);
else
return x;
}
}
}
}
}
}
strong_alias (__ieee754_remainder, __remainder_finite)

69
sysdeps/ieee754/dbl-64/e_sinh.c
View File

@ -38,43 +38,50 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
static const double one = 1.0, shuge = 1.0e307;
double
__ieee754_sinh(double x)
__ieee754_sinh (double x)
{
double t,w,h;
int32_t ix,jx;
u_int32_t lx;
double t, w, h;
int32_t ix, jx;
u_int32_t lx;
/* High word of |x|. */
GET_HIGH_WORD(jx,x);
ix = jx&0x7fffffff;
/* High word of |x|. */
GET_HIGH_WORD (jx, x);
ix = jx & 0x7fffffff;
/* x is INF or NaN */
if(__builtin_expect(ix>=0x7ff00000, 0)) return x+x;
/* x is INF or NaN */
if (__builtin_expect (ix >= 0x7ff00000, 0))
return x + x;
h = 0.5;
if (jx<0) h = -h;
/* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x40360000) { /* |x|<22 */
if (__builtin_expect(ix<0x3e300000, 0)) /* |x|<2**-28 */
if(shuge+x>one)
return x;/* sinh(tiny) = tiny with inexact */
t = __expm1(fabs(x));
if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
return h*(t+t/(t+one));
}
h = 0.5;
if (jx < 0)
h = -h;
/* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x40360000) /* |x|<22 */
{
if (__builtin_expect (ix < 0x3e300000, 0)) /* |x|<2**-28 */
if (shuge + x > one)
return x;
/* sinh(tiny) = tiny with inexact */
t = __expm1 (fabs (x));
if (ix < 0x3ff00000)
return h * (2.0 * t - t * t / (t + one));
return h * (t + t / (t + one));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862e42) return h*__ieee754_exp(fabs(x));
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862e42)
return h * __ieee754_exp (fabs (x));
/* |x| in [log(maxdouble), overflowthresold] */
GET_LOW_WORD(lx,x);
if (ix<0x408633ce || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
w = __ieee754_exp(0.5*fabs(x));
t = h*w;
return t*w;
}
/* |x| in [log(maxdouble), overflowthresold] */
GET_LOW_WORD (lx, x);
if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (u_int32_t) 0x8fb9f87d)))
{
w = __ieee754_exp (0.5 * fabs (x));
t = h * w;
return t * w;
}
/* |x| > overflowthresold, sinh(x) overflow */
return x*shuge;
/* |x| > overflowthresold, sinh(x) overflow */
return x * shuge;
}
strong_alias (__ieee754_sinh, __sinh_finite)

95
sysdeps/ieee754/dbl-64/e_sqrt.c
View File

@ -44,58 +44,67 @@
/* it computes the correctly rounded (to nearest) value of square */
/* root of x. */
/*********************************************************************/
double __ieee754_sqrt(double x) {
double
__ieee754_sqrt (double x)
{
#include "uroot.h"
static const double
rt0 = 9.99999999859990725855365213134618E-01,
rt1 = 4.99999999495955425917856814202739E-01,
rt2 = 3.75017500867345182581453026130850E-01,
rt3 = 3.12523626554518656309172508769531E-01;
static const double big = 134217728.0;
double y,t,del,res,res1,hy,z,zz,p,hx,tx,ty,s;
mynumber a,c={{0,0}};
static const double big = 134217728.0;
double y, t, del, res, res1, hy, z, zz, p, hx, tx, ty, s;
mynumber a, c = { { 0, 0 } };
int4 k;
a.x=x;
k=a.i[HIGH_HALF];
a.i[HIGH_HALF]=(k&0x001fffff)|0x3fe00000;
t=inroot[(k&0x001fffff)>>14];
s=a.x;
a.x = x;
k = a.i[HIGH_HALF];
a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000;
t = inroot[(k & 0x001fffff) >> 14];
s = a.x;
/*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
if (k>0x000fffff && k<0x7ff00000) {
fenv_t env;
libc_feholdexcept (&env);
double ret;
y=1.0-t*(t*s);
t=t*(rt0+y*(rt1+y*(rt2+y*rt3)));
c.i[HIGH_HALF]=0x20000000+((k&0x7fe00000)>>1);
y=t*s;
hy=(y+big)-big;
del=0.5*t*((s-hy*hy)-(y-hy)*(y+hy));
res=y+del;
if (res == (res+1.002*((y-res)+del))) ret = res*c.x;
else {
res1=res+1.5*((y-res)+del);
EMULV(res,res1,z,zz,p,hx,tx,hy,ty); /* (z+zz)=res*res1 */
ret = ((((z-s)+zz)<0)?max(res,res1):min(res,res1))*c.x;
if (k > 0x000fffff && k < 0x7ff00000)
{
fenv_t env;
libc_feholdexcept (&env);
double ret;
y = 1.0 - t * (t * s);
t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
c.i[HIGH_HALF] = 0x20000000 + ((k & 0x7fe00000) >> 1);
y = t * s;
hy = (y + big) - big;
del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy));
res = y + del;
if (res == (res + 1.002 * ((y - res) + del)))
ret = res * c.x;
else
{
res1 = res + 1.5 * ((y - res) + del);
EMULV (res, res1, z, zz, p, hx, tx, hy, ty); /* (z+zz)=res*res1 */
ret = ((((z - s) + zz) < 0) ? max (res, res1) :
min (res, res1)) * c.x;
}
math_force_eval (ret);
libc_fesetenv (&env);
if (x / ret != ret)
{
double force_inexact = 1.0 / 3.0;
math_force_eval (force_inexact);
}
/* Otherwise (x / ret == ret), either the square root was exact or
the division was inexact. */
return ret;
}
else
{
if ((k & 0x7ff00000) == 0x7ff00000)
return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
if (x == 0)
return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */
if (k < 0)
return (x - x) / (x - x); /* sqrt(-ve)=sNaN */
return tm256.x * __ieee754_sqrt (x * t512.x);
}
math_force_eval (ret);
libc_fesetenv (&env);
if (x / ret != ret)
{
double force_inexact = 1.0 / 3.0;
math_force_eval (force_inexact);
}
/* Otherwise (x / ret == ret), either the square root was exact or
the division was inexact. */
return ret;
}
else {
if ((k & 0x7ff00000) == 0x7ff00000)
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
if (x==0) return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */
if (k<0) return (x-x)/(x-x); /* sqrt(-ve)=sNaN */
return tm256.x*__ieee754_sqrt(x*t512.x);
}
}
strong_alias (__ieee754_sqrt, __sqrt_finite)

129
sysdeps/ieee754/dbl-64/halfulp.c
View File

@ -44,86 +44,109 @@
#endif
static const int4 tab54[32] = {
262143, 11585, 1782, 511, 210, 107, 63, 42,
30, 22, 17, 14, 12, 10, 9, 7,
7, 6, 5, 5, 5, 4, 4, 4,
3, 3, 3, 3, 3, 3, 3, 3 };
262143, 11585, 1782, 511, 210, 107, 63, 42,
30, 22, 17, 14, 12, 10, 9, 7,
7, 6, 5, 5, 5, 4, 4, 4,
3, 3, 3, 3, 3, 3, 3, 3
};
double
SECTION
__halfulp(double x, double y)
__halfulp (double x, double y)
{
mynumber v;
double z,u,uu;
double z, u, uu;
#ifndef DLA_FMS
double j1,j2,j3,j4,j5;
double j1, j2, j3, j4, j5;
#endif
int4 k,l,m,n;
if (y <= 0) { /*if power is negative or zero */
v.x = y;
if (v.i[LOW_HALF] != 0) return -10.0;
v.x = x;
if (v.i[LOW_HALF] != 0) return -10.0;
if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */
z = (double) k;
return (z*y == -1075.0)?0: -10.0;
}
/* if y > 0 */
int4 k, l, m, n;
if (y <= 0) /*if power is negative or zero */
{
v.x = y;
if (v.i[LOW_HALF] != 0)
return -10.0;
v.x = x;
if (v.i[LOW_HALF] != 0)
return -10.0;
if ((v.i[HIGH_HALF] & 0x000fffff) != 0)
return -10; /* if x =2 ^ n */
k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */
z = (double) k;
return (z * y == -1075.0) ? 0 : -10.0;
}
/* if y > 0 */
v.x = y;
if (v.i[LOW_HALF] != 0) return -10.0;
if (v.i[LOW_HALF] != 0)
return -10.0;
v.x=x;
/* case where x = 2**n for some integer n */
if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
k=(v.i[HIGH_HALF]>>20)-1023;
return (((double) k)*y == -1075.0)?0:-10.0;
}
v.x = x;
/* case where x = 2**n for some integer n */
if (((v.i[HIGH_HALF] & 0x000fffff) | v.i[LOW_HALF]) == 0)
{
k = (v.i[HIGH_HALF] >> 20) - 1023;
return (((double) k) * y == -1075.0) ? 0 : -10.0;
}
v.x = y;
k = v.i[HIGH_HALF];
m = k<<12;
m = k << 12;
l = 0;
while (m)
{m = m<<1; l++; }
n = (k&0x000fffff)|0x00100000;
n = n>>(20-l); /* n is the odd integer of y */
k = ((k>>20) -1023)-l; /* y = n*2**k */
if (k>5) return -10.0;
if (k>0) for (;k>0;k--) n *= 2;
if (n > 34) return -10.0;
{
m = m << 1; l++;
}
n = (k & 0x000fffff) | 0x00100000;
n = n >> (20 - l); /* n is the odd integer of y */
k = ((k >> 20) - 1023) - l; /* y = n*2**k */
if (k > 5)
return -10.0;
if (k > 0)
for (; k > 0; k--)
n *= 2;
if (n > 34)
return -10.0;
k = -k;
if (k>5) return -10.0;
if (k > 5)
return -10.0;
/* now treat x */
while (k>0) {
z = __ieee754_sqrt(x);
EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
if (((u-x)+uu) != 0) break;
x = z;
k--;
}
if (k) return -10.0;
/* now treat x */
while (k > 0)
{
z = __ieee754_sqrt (x);
EMULV (z, z, u, uu, j1, j2, j3, j4, j5);
if (((u - x) + uu) != 0)
break;
x = z;
k--;
}
if (k)
return -10.0;
/* it is impossible that n == 2, so the mantissa of x must be short */
v.x = x;
if (v.i[LOW_HALF]) return -10.0;
if (v.i[LOW_HALF])
return -10.0;
k = v.i[HIGH_HALF];
m = k<<12;
m = k << 12;
l = 0;
while (m) {m = m<<1; l++; }
m = (k&0x000fffff)|0x00100000;
m = m>>(20-l); /* m is the odd integer of x */
while (m)
{
m = m << 1; l++;
}
m = (k & 0x000fffff) | 0x00100000;
m = m >> (20 - l); /* m is the odd integer of x */
/* now check whether the length of m**n is at most 54 bits */
/* now check whether the length of m**n is at most 54 bits */
if (m > tab54[n-3]) return -10.0;
if (m > tab54[n - 3])
return -10.0;
/* yes, it is - now compute x**n by simple multiplications */
/* yes, it is - now compute x**n by simple multiplications */
u = x;
for (k=1;k<n;k++) u = u*x;
for (k = 1; k < n; k++)
u = u * x;
return u;
}

343
sysdeps/ieee754/dbl-64/k_rem_pio2.c
View File

@ -147,166 +147,215 @@ static const double PIo2[] = {
};
static const double
zero = 0.0,
one = 1.0,
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
zero = 0.0,
one = 1.0,
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
int
__kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
const int32_t *ipio2)
{
int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
double z,fw,f[20],fq[20],q[20];
int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
double z, fw, f[20], fq[20], q[20];
/* initialize jk*/
jk = init_jk[prec];
jp = jk;
/* initialize jk*/
jk = init_jk[prec];
jp = jk;
/* determine jx,jv,q0, note that 3>q0 */
jx = nx-1;
jv = (e0-3)/24; if(jv<0) jv=0;
q0 = e0-24*(jv+1);
/* determine jx,jv,q0, note that 3>q0 */
jx = nx - 1;
jv = (e0 - 3) / 24; if (jv < 0)
jv = 0;
q0 = e0 - 24 * (jv + 1);
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv-jx; m = jx+jk;
for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv - jx; m = jx + jk;
for (i = 0; i <= m; i++, j++)
f[i] = (j < 0) ? zero : (double) ipio2[j];
/* compute q[0],q[1],...q[jk] */
for (i=0;i<=jk;i++) {
for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
}
/* compute q[0],q[1],...q[jk] */
for (i = 0; i <= jk; i++)
{
for (j = 0, fw = 0.0; j <= jx; j++)
fw += x[j] * f[jx + i - j];
q[i] = fw;
}
jz = jk;
jz = jk;
recompute:
/* distill q[] into iq[] reversingly */
for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
fw = (double)((int32_t)(twon24* z));
iq[i] = (int32_t)(z-two24*fw);
z = q[j-1]+fw;
}
/* distill q[] into iq[] reversingly */
for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
{
fw = (double) ((int32_t) (twon24 * z));
iq[i] = (int32_t) (z - two24 * fw);
z = q[j - 1] + fw;
}
/* compute n */
z = __scalbn(z,q0); /* actual value of z */
z -= 8.0*__floor(z*0.125); /* trim off integer >= 8 */
n = (int32_t) z;
z -= (double)n;
ih = 0;
if(q0>0) { /* need iq[jz-1] to determine n */
i = (iq[jz-1]>>(24-q0)); n += i;
iq[jz-1] -= i<<(24-q0);
ih = iq[jz-1]>>(23-q0);
}
else if(q0==0) ih = iq[jz-1]>>23;
else if(z>=0.5) ih=2;
/* compute n */
z = __scalbn (z, q0); /* actual value of z */
z -= 8.0 * __floor (z * 0.125); /* trim off integer >= 8 */
n = (int32_t) z;
z -= (double) n;
ih = 0;
if (q0 > 0) /* need iq[jz-1] to determine n */
{
i = (iq[jz - 1] >> (24 - q0)); n += i;
iq[jz - 1] -= i << (24 - q0);
ih = iq[jz - 1] >> (23 - q0);
}
else if (q0 == 0)
ih = iq[jz - 1] >> 23;
else if (z >= 0.5)
ih = 2;
if(ih>0) { /* q > 0.5 */
n += 1; carry = 0;
for(i=0;i<jz ;i++) { /* compute 1-q */
j = iq[i];
if(carry==0) {
if(j!=0) {
carry = 1; iq[i] = 0x1000000- j;
}
} else iq[i] = 0xffffff - j;
}
if(q0>0) { /* rare case: chance is 1 in 12 */
switch(q0) {
case 1:
iq[jz-1] &= 0x7fffff; break;
case 2:
iq[jz-1] &= 0x3fffff; break;
}
}
if(ih==2) {
z = one - z;
if(carry!=0) z -= __scalbn(one,q0);
}
}
/* check if recomputation is needed */
if(z==zero) {
j = 0;
for (i=jz-1;i>=jk;i--) j |= iq[i];
if(j==0) { /* need recomputation */
for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
f[jx+i] = (double) ipio2[jv+i];
for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
q[i] = fw;
if (ih > 0) /* q > 0.5 */
{
n += 1; carry = 0;
for (i = 0; i < jz; i++) /* compute 1-q */
{
j = iq[i];
if (carry == 0)
{
if (j != 0)
{
carry = 1; iq[i] = 0x1000000 - j;
}
jz += k;
goto recompute;
}
else
iq[i] = 0xffffff - j;
}
/* chop off zero terms */
if(z==0.0) {
jz -= 1; q0 -= 24;
while(iq[jz]==0) { jz--; q0-=24;}
} else { /* break z into 24-bit if necessary */
z = __scalbn(z,-q0);
if(z>=two24) {
fw = (double)((int32_t)(twon24*z));
iq[jz] = (int32_t)(z-two24*fw);
jz += 1; q0 += 24;
iq[jz] = (int32_t) fw;
} else iq[jz] = (int32_t) z ;
}
/* convert integer "bit" chunk to floating-point value */
fw = __scalbn(one,q0);
for(i=jz;i>=0;i--) {
q[i] = fw*(double)iq[i]; fw*=twon24;
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */
for(i=jz;i>=0;i--) {
for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
fq[jz-i] = fw;
}
/* compress fq[] into y[] */
switch(prec) {
case 0:
fw = 0.0;
for (i=jz;i>=0;i--) fw += fq[i];
y[0] = (ih==0)? fw: -fw;
break;
if (q0 > 0) /* rare case: chance is 1 in 12 */
{
switch (q0)
{
case 1:
case 2:;
#if __FLT_EVAL_METHOD__ != 0
volatile
#endif
double fv = 0.0;
for (i=jz;i>=0;i--) fv += fq[i];
y[0] = (ih==0)? fv: -fv;
fv = fq[0]-fv;
for (i=1;i<=jz;i++) fv += fq[i];
y[1] = (ih==0)? fv: -fv;
break;
case 3: /* painful */
for (i=jz;i>0;i--) {
#if __FLT_EVAL_METHOD__ != 0
volatile
#endif
double fv = (double)(fq[i-1]+fq[i]);
fq[i] += fq[i-1]-fv;
fq[i-1] = fv;
}
for (i=jz;i>1;i--) {
#if __FLT_EVAL_METHOD__ != 0
volatile
#endif
double fv = (double)(fq[i-1]+fq[i]);
fq[i] += fq[i-1]-fv;
fq[i-1] = fv;
}
for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
if(ih==0) {
y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
} else {
y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
}
iq[jz - 1] &= 0x7fffff; break;
case 2:
iq[jz - 1] &= 0x3fffff; break;
}
}
return n&7;
if (ih == 2)
{
z = one - z;
if (carry != 0)
z -= __scalbn (one, q0);
}
}
/* check if recomputation is needed */
if (z == zero)
{
j = 0;
for (i = jz - 1; i >= jk; i--)
j |= iq[i];
if (j == 0) /* need recomputation */
{
for (k = 1; iq[jk - k] == 0; k++)
; /* k = no. of terms needed */
for (i = jz + 1; i <= jz + k; i++) /* add q[jz+1] to q[jz+k] */
{
f[jx + i] = (double) ipio2[jv + i];
for (j = 0, fw = 0.0; j <= jx; j++)
fw += x[j] * f[jx + i - j];
q[i] = fw;
}
jz += k;
goto recompute;
}
}
/* chop off zero terms */
if (z == 0.0)
{
jz -= 1; q0 -= 24;
while (iq[jz] == 0)
{
jz--; q0 -= 24;
}
}
else /* break z into 24-bit if necessary */
{
z = __scalbn (z, -q0);
if (z >= two24)
{
fw = (double) ((int32_t) (twon24 * z));
iq[jz] = (int32_t) (z - two24 * fw);
jz += 1; q0 += 24;
iq[jz] = (int32_t) fw;
}
else
iq[jz] = (int32_t) z;
}
/* convert integer "bit" chunk to floating-point value */
fw = __scalbn (one, q0);
for (i = jz; i >= 0; i--)
{
q[i] = fw * (double) iq[i]; fw *= twon24;
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */
for (i = jz; i >= 0; i--)
{
for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
fw += PIo2[k] * q[i + k];
fq[jz - i] = fw;
}
/* compress fq[] into y[] */
switch (prec)
{
case 0:
fw = 0.0;
for (i = jz; i >= 0; i--)
fw += fq[i];
y[0] = (ih == 0) ? fw : -fw;
break;
case 1:
case 2:;
#if __FLT_EVAL_METHOD__ != 0
volatile
#endif
double fv = 0.0;
for (i = jz; i >= 0; i--)
fv += fq[i];
y[0] = (ih == 0) ? fv : -fv;
fv = fq[0] - fv;
for (i = 1; i <= jz; i++)
fv += fq[i];
y[1] = (ih == 0) ? fv : -fv;
break;
case 3: /* painful */
for (i = jz; i > 0; i--)
{
#if __FLT_EVAL_METHOD__ != 0
volatile
#endif
double fv = (double) (fq[i - 1] + fq[i]);
fq[i] += fq[i - 1] - fv;
fq[i - 1] = fv;
}
for (i = jz; i > 1; i--)
{
#if __FLT_EVAL_METHOD__ != 0
volatile
#endif
double fv = (double) (fq[i - 1] + fq[i]);
fq[i] += fq[i - 1] - fv;
fq[i - 1] = fv;
}
for (fw = 0.0, i = jz; i >= 2; i--)
fw += fq[i];
if (ih == 0)
{
y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
}
else
{
y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
}
}
return n & 7;
}

20
sysdeps/ieee754/dbl-64/mpa-arch.h
View File

@ -22,15 +22,15 @@ typedef int64_t mantissa_store_t;
#define TWOPOW(i) (1L << i)
#define RADIX_EXP 24
#define RADIX TWOPOW (RADIX_EXP) /* 2^24 */
#define RADIX TWOPOW (RADIX_EXP) /* 2^24 */
/* Divide D by RADIX and put the remainder in R. D must be a non-negative
integral value. */
#define DIV_RADIX(d, r) \
({ \
r = d & (RADIX - 1); \
d >>= RADIX_EXP; \
})
({ \
r = d & (RADIX - 1); \
d >>= RADIX_EXP; \
})
/* Put the integer component of a double X in R and retain the fraction in
X. This is used in extracting mantissa digits for MP_NO by using the
@ -38,10 +38,10 @@ typedef int64_t mantissa_store_t;
digit and then scaling by RADIX to get the next mantissa digit in the same
manner. */
#define INTEGER_OF(x, i) \
({ \
i = (mantissa_t) x; \
x -= i; \
})
({ \
i = (mantissa_t) x; \
x -= i; \
})
/* Align IN down to F. The code assumes that F is a power of two. */
#define ALIGN_DOWN_TO(in, f) ((in) & -(f))
#define ALIGN_DOWN_TO(in, f) ((in) & - (f))

29
sysdeps/ieee754/dbl-64/mpa.c
View File

@ -50,8 +50,8 @@
#endif
#ifndef NO__CONST
const mp_no mpone = {1, {1.0, 1.0}};
const mp_no mptwo = {1, {1.0, 2.0}};
const mp_no mpone = { 1, { 1.0, 1.0 } };
const mp_no mptwo = { 1, { 1.0, 2.0 } };
#endif
#ifndef NO___ACR
@ -123,7 +123,7 @@ __cpy (const mp_no *x, mp_no *y, int p)
static void
norm (const mp_no *x, double *y, int p)
{
#define R RADIXI
# define R RADIXI
long i;
double c;
mantissa_t a, u, v, z[5];
@ -140,7 +140,7 @@ norm (const mp_no *x, double *y, int p)
}
else
{
for (a = 1, z[1] = X[1]; z[1] < TWO23;)
for (a = 1, z[1] = X[1]; z[1] < TWO23; )
{
a *= 2;
z[1] *= 2;
@ -188,7 +188,7 @@ norm (const mp_no *x, double *y, int p)
c *= RADIXI;
*y = c;
#undef R
# undef R
}
/* Convert a multiple precision number *X into a double precision
@ -201,7 +201,7 @@ denorm (const mp_no *x, double *y, int p)
double c;
mantissa_t u, z[5];
#define R RADIXI
# define R RADIXI
if (EX < -44 || (EX == -44 && X[1] < TWO5))
{
*y = 0;
@ -298,7 +298,7 @@ denorm (const mp_no *x, double *y, int p)
c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
*y = c * TWOM1032;
#undef R
# undef R
}
/* Convert multiple precision number *X into double precision number *Y. The
@ -394,7 +394,7 @@ add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
zk = 1;
}
else
{
{
Z[k--] = zk;
zk = 0;
}
@ -409,7 +409,7 @@ add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
zk = 1;
}
else
{
{
Z[k--] = zk;
zk = 0;
}
@ -471,7 +471,7 @@ sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
zk = -1;
}
else
{
{
Z[k--] = zk;
zk = 0;
}
@ -487,18 +487,19 @@ sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
zk = -1;
}
else
{
{
Z[k--] = zk;
zk = 0;
}
}
/* Normalize. */
for (i = 1; Z[i] == 0; i++);
for (i = 1; Z[i] == 0; i++)
;
EZ = EZ - i + 1;
for (k = 1; i <= p2 + 1;)
for (k = 1; i <= p2 + 1; )
Z[k++] = Z[i++];
for (; k <= p2;)
for (; k <= p2; )
Z[k++] = 0;
}

1
sysdeps/ieee754/dbl-64/mpatan.c
View File

@ -43,7 +43,6 @@ void
SECTION
__mpatan (mp_no *x, mp_no *y, int p)
{
int i, m, n;
double dx;
mp_no mptwoim1 =

2
sysdeps/ieee754/dbl-64/mpn2dbl.c
View File

@ -40,7 +40,7 @@ __mpn_construct_double (mp_srcptr frac_ptr, int expt, int negative)
u.ieee.mantissa0 = (frac_ptr[0] >> 32) & (((mp_limb_t) 1
<< (DBL_MANT_DIG - 32)) - 1);
#else
#error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
# error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
#endif
return u.d;

1
sysdeps/ieee754/dbl-64/mptan.c
View File

@ -44,7 +44,6 @@ void
SECTION
__mptan (double x, mp_no *mpy, int p)
{
int n;
mp_no mpw, mpc, mps;

9
sysdeps/ieee754/dbl-64/mydefs.h
View File

@ -28,10 +28,9 @@
#define MY_H
typedef int int4;
typedef union {int4 i[2]; double x;} mynumber;
#define ABS(x) (((x)>0)?(x):-(x))
#define max(x,y) (((y)>(x))?(y):(x))
#define min(x,y) (((y)<(x))?(y):(x))
typedef union { int4 i[2]; double x; } mynumber;
#define ABS(x) (((x) > 0) ? (x) : -(x))
#define max(x, y) (((y) > (x)) ? (y) : (x))
#define min(x, y) (((y) < (x)) ? (y) : (x))
#endif

54
sysdeps/ieee754/dbl-64/s_asinh.c
View File

@ -25,33 +25,43 @@
#include <math_private.h>
static const double
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300;
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge = 1.00000000000000000000e+300;
double
__asinh(double x)
__asinh (double x)
{
double w;
int32_t hx,ix;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
if(__builtin_expect(ix< 0x3e300000, 0)) { /* |x|<2**-28 */
if(huge+x>one) return x; /* return x inexact except 0 */
double w;
int32_t hx, ix;
GET_HIGH_WORD (hx, x);
ix = hx & 0x7fffffff;
if (__builtin_expect (ix < 0x3e300000, 0)) /* |x|<2**-28 */
{
if (huge + x > one)
return x; /* return x inexact except 0 */
}
if (__builtin_expect (ix > 0x41b00000, 0)) /* |x| > 2**28 */
{
if (ix >= 0x7ff00000)
return x + x; /* x is inf or NaN */
w = __ieee754_log (fabs (x)) + ln2;
}
else
{
double xa = fabs (x);
if (ix > 0x40000000) /* 2**28 > |x| > 2.0 */
{
w = __ieee754_log (2.0 * xa + one / (__ieee754_sqrt (xa * xa + one) +
xa));
}
if(__builtin_expect(ix>0x41b00000, 0)) { /* |x| > 2**28 */
if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */
w = __ieee754_log(fabs(x))+ln2;
} else {
double xa = fabs(x);
if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
w = __ieee754_log(2.0*xa+one/(__ieee754_sqrt(xa*xa+one)+xa));
} else { /* 2.0 > |x| > 2**-28 */
double t = xa*xa;
w =__log1p(xa+t/(one+__ieee754_sqrt(one+t)));
}
else /* 2.0 > |x| > 2**-28 */
{
double t = xa * xa;
w = __log1p (xa + t / (one + __ieee754_sqrt (one + t)));
}
return __copysign(w, x);
}
return __copysign (w, x);
}
weak_alias (__asinh, asinh)
#ifdef NO_LONG_DOUBLE

14
sysdeps/ieee754/dbl-64/s_atan.c
View File

@ -61,7 +61,7 @@ double
atan (double x)
{
double cor, s1, ss1, s2, ss2, t1, t2, t3, t7, t8, t9, t10, u, u2, u3,
v, vv, w, ww, y, yy, z, zz;
v, vv, w, ww, y, yy, z, zz;
#ifndef DLA_FMS
double t4, t5, t6;
#endif
@ -191,17 +191,17 @@ atan (double x)
yy = cij[i][4].d + z * yy;
yy = cij[i][3].d + z * yy;
yy = cij[i][2].d + z * yy;
yy = HPI1 - z * yy;
yy = HPI1 - z * yy;
t1 = HPI - cij[i][1].d;
if (i < 112)
u3 = U31; /* w < 1/2 */
u3 = U31; /* w < 1/2 */
else
u3 = U32; /* w >= 1/2 */
u3 = U32; /* w >= 1/2 */
if ((y = t1 + (yy - u3)) == t1 + (yy + u3))
return __signArctan (x, y);
DIV2 (1 , 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8, t9,
DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
t1 = w - hij[i][0].d;
EADD (t1, ww, z, zz);
@ -230,7 +230,7 @@ atan (double x)
else
{
if (u < E)
{ /* D <= u < E */
{ /* D <= u < E */
w = 1 / u;
v = w * w;
EMULV (w, u, t1, t2, t3, t4, t5, t6, t7);
@ -248,7 +248,7 @@ atan (double x)
if ((y = t3 + (yy - U4)) == t3 + (yy + U4))
return __signArctan (x, y);
DIV2 (1 , 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8,
DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8,
t9, t10);
MUL2 (w, ww, w, ww, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);

19
sysdeps/ieee754/dbl-64/s_cbrt.c
View File

@ -51,16 +51,17 @@ __cbrt (double x)
if (xe == 0 && fpclassify (x) <= FP_ZERO)
return x + x;
u = (0.354895765043919860
+ ((1.50819193781584896
+ ((-2.11499494167371287
+ ((2.44693122563534430
+ ((-1.83469277483613086
+ (0.784932344976639262 - 0.145263899385486377 * xm) * xm)
u = (0.354895765043919860
+ ((1.50819193781584896
+ ((-2.11499494167371287
+ ((2.44693122563534430
+ ((-1.83469277483613086
+ (0.784932344976639262 - 0.145263899385486377 * xm)
* xm)
* xm))
* xm))
* xm))
* xm))
* xm));
* xm))
* xm));
t2 = u * u * u;

91
sysdeps/ieee754/dbl-64/s_ceil.c
View File

@ -25,44 +25,67 @@
static const double huge = 1.0e300;
double
__ceil(double x)
__ceil (double x)
{
int32_t i0,i1,j0;
u_int32_t i,j;
EXTRACT_WORDS(i0,i1,x);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) { /* raise inexact if x != 0 */
math_force_eval(huge+x);
/* return 0*sign(x) if |x|<1 */
if(i0<0) {i0=0x80000000;i1=0;}
else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) return x; /* x is integral */
math_force_eval(huge+x); /* raise inexact flag */
if(i0>0) i0 += (0x00100000)>>j0;
i0 &= (~i); i1=0;
int32_t i0, i1, j0;
u_int32_t i, j;
EXTRACT_WORDS (i0, i1, x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20)
{
if (j0 < 0) /* raise inexact if x != 0 */
{
math_force_eval (huge + x);
/* return 0*sign(x) if |x|<1 */
if (i0 < 0)
{
i0 = 0x80000000; i1 = 0;
}
} else if (j0>51) {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) return x; /* x is integral */
math_force_eval(huge+x); /* raise inexact flag */
if(i0>0) {
if(j0==20) i0+=1;
else {
j = i1 + (1<<(52-j0));
if(j<i1) i0+=1; /* got a carry */
i1 = j;
}
else if ((i0 | i1) != 0)
{
i0 = 0x3ff00000; i1 = 0;
}
i1 &= (~i);
}
INSERT_WORDS(x,i0,i1);
return x;
else
{
i = (0x000fffff) >> j0;
if (((i0 & i) | i1) == 0)
return x; /* x is integral */
math_force_eval (huge + x); /* raise inexact flag */
if (i0 > 0)
i0 += (0x00100000) >> j0;
i0 &= (~i); i1 = 0;
}
}
else if (j0 > 51)
{
if (j0 == 0x400)
return x + x; /* inf or NaN */
else
return x; /* x is integral */
}
else
{
i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0)
return x; /* x is integral */
math_force_eval (huge + x); /* raise inexact flag */
if (i0 > 0)
{
if (j0 == 20)
i0 += 1;
else
{
j = i1 + (1 << (52 - j0));
if (j < i1)
i0 += 1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
INSERT_WORDS (x, i0, i1);
return x;
}
#ifndef __ceil
weak_alias (__ceil, ceil)

13
sysdeps/ieee754/dbl-64/s_copysign.c
View File

@ -23,13 +23,14 @@ static char rcsid[] = "$NetBSD: s_copysign.c,v 1.8 1995/05/10 20:46:57 jtc Exp $
#include <math.h>
#include <math_private.h>
double __copysign(double x, double y)
double
__copysign (double x, double y)
{
u_int32_t hx,hy;
GET_HIGH_WORD(hx,x);
GET_HIGH_WORD(hy,y);
SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
return x;
u_int32_t hx, hy;
GET_HIGH_WORD (hx, x);
GET_HIGH_WORD (hy, y);
SET_HIGH_WORD (x, (hx & 0x7fffffff) | (hy & 0x80000000));
return x;
}
weak_alias (__copysign, copysign)
#ifdef NO_LONG_DOUBLE

475
sysdeps/ieee754/dbl-64/s_erf.c
View File

@ -116,157 +116,176 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
#include <math_private.h>
static const double
tiny = 1e-300,
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
tiny = 1e-300,
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp[] = {1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
-3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
-2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
-5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
-2.37630166566501626084e-05}, /* 0xBEF8EAD6, 0x120016AC */
qq[] = {0.0, 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
-3.96022827877536812320e-06}, /* 0xBED09C43, 0x42A26120 */
efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp[] = { 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
-3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
-2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
-5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
-2.37630166566501626084e-05 }, /* 0xBEF8EAD6, 0x120016AC */
qq[] = { 0.0, 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
-3.96022827877536812320e-06 }, /* 0xBED09C43, 0x42A26120 */
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
pa[] = {-2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
-3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
-1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
-2.16637559486879084300e-03}, /* 0xBF61BF38, 0x0A96073F */
qa[] = {0.0, 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
1.19844998467991074170e-02}, /* 0x3F888B54, 0x5735151D */
pa[] = { -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
-3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
-1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
-2.16637559486879084300e-03 }, /* 0xBF61BF38, 0x0A96073F */
qa[] = { 0.0, 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
1.19844998467991074170e-02 }, /* 0x3F888B54, 0x5735151D */
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
ra[] = {-9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
-6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
-1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
-6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
-1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
-1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
-8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
-9.81432934416914548592e+00}, /* 0xC023A0EF, 0xC69AC25C */
sa[] = {0.0,1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
-6.04244152148580987438e-02}, /* 0xBFAEEFF2, 0xEE749A62 */
ra[] = { -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
-6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
-1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
-6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
-1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
-1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
-8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
-9.81432934416914548592e+00 }, /* 0xC023A0EF, 0xC69AC25C */
sa[] = { 0.0, 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
-6.04244152148580987438e-02 }, /* 0xBFAEEFF2, 0xEE749A62 */
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
rb[] = {-9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
-7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
-1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
-1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
-6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
-1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
-4.83519191608651397019e+02}, /* 0xC07E384E, 0x9BDC383F */
sb[] = {0.0,3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
-2.24409524465858183362e+01}; /* 0xC03670E2, 0x42712D62 */
rb[] = { -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
-7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
-1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
-1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
-6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
-1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
-4.83519191608651397019e+02 }, /* 0xC07E384E, 0x9BDC383F */
sb[] = { 0.0, 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
-2.24409524465858183362e+01 }; /* 0xC03670E2, 0x42712D62 */
double __erf(double x)
double
__erf (double x)
{
int32_t hx,ix,i;
double R,S,P,Q,s,y,z,r;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erf(nan)=nan */
i = ((u_int32_t)hx>>31)<<1;
return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
}
int32_t hx, ix, i;
double R, S, P, Q, s, y, z, r;
GET_HIGH_WORD (hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) /* erf(nan)=nan */
{
i = ((u_int32_t) hx >> 31) << 1;
return (double) (1 - i) + one / x; /* erf(+-inf)=+-1 */
}
if(ix < 0x3feb0000) { /* |x|<0.84375 */
double r1,r2,s1,s2,s3,z2,z4;
if(ix < 0x3e300000) { /* |x|<2**-28 */
if (ix < 0x00800000)
return 0.125*(8.0*x+efx8*x); /*avoid underflow */
return x + efx*x;
}
z = x*x;
r1 = pp[0]+z*pp[1]; z2=z*z;
r2 = pp[2]+z*pp[3]; z4=z2*z2;
s1 = one+z*qq[1];
s2 = qq[2]+z*qq[3];
s3 = qq[4]+z*qq[5];
r = r1 + z2*r2 + z4*pp[4];
s = s1 + z2*s2 + z4*s3;
y = r/s;
return x + x*y;
if (ix < 0x3feb0000) /* |x|<0.84375 */
{
double r1, r2, s1, s2, s3, z2, z4;
if (ix < 0x3e300000) /* |x|<2**-28 */
{
if (ix < 0x00800000)
return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */
return x + efx * x;
}
if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
double s2,s4,s6,P1,P2,P3,P4,Q1,Q2,Q3,Q4;
s = fabs(x)-one;
P1 = pa[0]+s*pa[1]; s2=s*s;
Q1 = one+s*qa[1]; s4=s2*s2;
P2 = pa[2]+s*pa[3]; s6=s4*s2;
Q2 = qa[2]+s*qa[3];
P3 = pa[4]+s*pa[5];
Q3 = qa[4]+s*qa[5];
P4 = pa[6];
Q4 = qa[6];
P = P1 + s2*P2 + s4*P3 + s6*P4;
Q = Q1 + s2*Q2 + s4*Q3 + s6*Q4;
if(hx>=0) return erx + P/Q; else return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf>|x|>=6 */
if(hx>=0) return one-tiny; else return tiny-one;
}
x = fabs(x);
s = one/(x*x);
if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */
double R1,R2,R3,R4,S1,S2,S3,S4,s2,s4,s6,s8;
R1 = ra[0]+s*ra[1];s2 = s*s;
S1 = one+s*sa[1]; s4 = s2*s2;
R2 = ra[2]+s*ra[3];s6 = s4*s2;
S2 = sa[2]+s*sa[3];s8 = s4*s4;
R3 = ra[4]+s*ra[5];
S3 = sa[4]+s*sa[5];
R4 = ra[6]+s*ra[7];
S4 = sa[6]+s*sa[7];
R = R1 + s2*R2 + s4*R3 + s6*R4;
S = S1 + s2*S2 + s4*S3 + s6*S4 + s8*sa[8];
} else { /* |x| >= 1/0.35 */
double R1,R2,R3,S1,S2,S3,S4,s2,s4,s6;
R1 = rb[0]+s*rb[1];s2 = s*s;
S1 = one+s*sb[1]; s4 = s2*s2;
R2 = rb[2]+s*rb[3];s6 = s4*s2;
S2 = sb[2]+s*sb[3];
R3 = rb[4]+s*rb[5];
S3 = sb[4]+s*sb[5];
S4 = sb[6]+s*sb[7];
R = R1 + s2*R2 + s4*R3 + s6*rb[6];
S = S1 + s2*S2 + s4*S3 + s6*S4;
}
z = x;
SET_LOW_WORD(z,0);
r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
if(hx>=0) return one-r/x; else return r/x-one;
z = x * x;
r1 = pp[0] + z * pp[1]; z2 = z * z;
r2 = pp[2] + z * pp[3]; z4 = z2 * z2;
s1 = one + z * qq[1];
s2 = qq[2] + z * qq[3];
s3 = qq[4] + z * qq[5];
r = r1 + z2 * r2 + z4 * pp[4];
s = s1 + z2 * s2 + z4 * s3;
y = r / s;
return x + x * y;
}
if (ix < 0x3ff40000) /* 0.84375 <= |x| < 1.25 */
{
double s2, s4, s6, P1, P2, P3, P4, Q1, Q2, Q3, Q4;
s = fabs (x) - one;
P1 = pa[0] + s * pa[1]; s2 = s * s;
Q1 = one + s * qa[1]; s4 = s2 * s2;
P2 = pa[2] + s * pa[3]; s6 = s4 * s2;
Q2 = qa[2] + s * qa[3];
P3 = pa[4] + s * pa[5];
Q3 = qa[4] + s * qa[5];
P4 = pa[6];
Q4 = qa[6];
P = P1 + s2 * P2 + s4 * P3 + s6 * P4;
Q = Q1 + s2 * Q2 + s4 * Q3 + s6 * Q4;
if (hx >= 0)
return erx + P / Q;
else
return -erx - P / Q;
}
if (ix >= 0x40180000) /* inf>|x|>=6 */
{
if (hx >= 0)
return one - tiny;
else
return tiny - one;
}
x = fabs (x);
s = one / (x * x);
if (ix < 0x4006DB6E) /* |x| < 1/0.35 */
{
double R1, R2, R3, R4, S1, S2, S3, S4, s2, s4, s6, s8;
R1 = ra[0] + s * ra[1]; s2 = s * s;
S1 = one + s * sa[1]; s4 = s2 * s2;
R2 = ra[2] + s * ra[3]; s6 = s4 * s2;
S2 = sa[2] + s * sa[3]; s8 = s4 * s4;
R3 = ra[4] + s * ra[5];
S3 = sa[4] + s * sa[5];
R4 = ra[6] + s * ra[7];
S4 = sa[6] + s * sa[7];
R = R1 + s2 * R2 + s4 * R3 + s6 * R4;
S = S1 + s2 * S2 + s4 * S3 + s6 * S4 + s8 * sa[8];
}
else /* |x| >= 1/0.35 */
{
double R1, R2, R3, S1, S2, S3, S4, s2, s4, s6;
R1 = rb[0] + s * rb[1]; s2 = s * s;
S1 = one + s * sb[1]; s4 = s2 * s2;
R2 = rb[2] + s * rb[3]; s6 = s4 * s2;
S2 = sb[2] + s * sb[3];
R3 = rb[4] + s * rb[5];
S3 = sb[4] + s * sb[5];
S4 = sb[6] + s * sb[7];
R = R1 + s2 * R2 + s4 * R3 + s6 * rb[6];
S = S1 + s2 * S2 + s4 * S3 + s6 * S4;
}
z = x;
SET_LOW_WORD (z, 0);
r = __ieee754_exp (-z * z - 0.5625) *
__ieee754_exp ((z - x) * (z + x) + R / S);
if (hx >= 0)
return one - r / x;
else
return r / x - one;
}
weak_alias (__erf, erf)
#ifdef NO_LONG_DOUBLE
@ -274,93 +293,115 @@ strong_alias (__erf, __erfl)
weak_alias (__erf, erfl)
#endif
double __erfc(double x)
double
__erfc (double x)
{
int32_t hx,ix;
double R,S,P,Q,s,y,z,r;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (double)(((u_int32_t)hx>>31)<<1)+one/x;
}
int32_t hx, ix;
double R, S, P, Q, s, y, z, r;
GET_HIGH_WORD (hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) /* erfc(nan)=nan */
{ /* erfc(+-inf)=0,2 */
return (double) (((u_int32_t) hx >> 31) << 1) + one / x;
}
if(ix < 0x3feb0000) { /* |x|<0.84375 */
double r1,r2,s1,s2,s3,z2,z4;
if(ix < 0x3c700000) /* |x|<2**-56 */
return one-x;
z = x*x;
r1 = pp[0]+z*pp[1]; z2=z*z;
r2 = pp[2]+z*pp[3]; z4=z2*z2;
s1 = one+z*qq[1];
s2 = qq[2]+z*qq[3];
s3 = qq[4]+z*qq[5];
r = r1 + z2*r2 + z4*pp[4];
s = s1 + z2*s2 + z4*s3;
y = r/s;
if(hx < 0x3fd00000) { /* x<1/4 */
return one-(x+x*y);
} else {
r = x*y;
r += (x-half);
return half - r ;
}
if (ix < 0x3feb0000) /* |x|<0.84375 */
{
double r1, r2, s1, s2, s3, z2, z4;
if (ix < 0x3c700000) /* |x|<2**-56 */
return one - x;
z = x * x;
r1 = pp[0] + z * pp[1]; z2 = z * z;
r2 = pp[2] + z * pp[3]; z4 = z2 * z2;
s1 = one + z * qq[1];
s2 = qq[2] + z * qq[3];
s3 = qq[4] + z * qq[5];
r = r1 + z2 * r2 + z4 * pp[4];
s = s1 + z2 * s2 + z4 * s3;
y = r / s;
if (hx < 0x3fd00000) /* x<1/4 */
{
return one - (x + x * y);
}
if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
double s2,s4,s6,P1,P2,P3,P4,Q1,Q2,Q3,Q4;
s = fabs(x)-one;
P1 = pa[0]+s*pa[1]; s2=s*s;
Q1 = one+s*qa[1]; s4=s2*s2;
P2 = pa[2]+s*pa[3]; s6=s4*s2;
Q2 = qa[2]+s*qa[3];
P3 = pa[4]+s*pa[5];
Q3 = qa[4]+s*qa[5];
P4 = pa[6];
Q4 = qa[6];
P = P1 + s2*P2 + s4*P3 + s6*P4;
Q = Q1 + s2*Q2 + s4*Q3 + s6*Q4;
if(hx>=0) {
z = one-erx; return z - P/Q;
} else {
z = erx+P/Q; return one+z;
}
else
{
r = x * y;
r += (x - half);
return half - r;
}
if (ix < 0x403c0000) { /* |x|<28 */
x = fabs(x);
s = one/(x*x);
if(ix< 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
double R1,R2,R3,R4,S1,S2,S3,S4,s2,s4,s6,s8;
R1 = ra[0]+s*ra[1];s2 = s*s;
S1 = one+s*sa[1]; s4 = s2*s2;
R2 = ra[2]+s*ra[3];s6 = s4*s2;
S2 = sa[2]+s*sa[3];s8 = s4*s4;
R3 = ra[4]+s*ra[5];
S3 = sa[4]+s*sa[5];
R4 = ra[6]+s*ra[7];
S4 = sa[6]+s*sa[7];
R = R1 + s2*R2 + s4*R3 + s6*R4;
S = S1 + s2*S2 + s4*S3 + s6*S4 + s8*sa[8];
} else { /* |x| >= 1/.35 ~ 2.857143 */
double R1,R2,R3,S1,S2,S3,S4,s2,s4,s6;
if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
R1 = rb[0]+s*rb[1];s2 = s*s;
S1 = one+s*sb[1]; s4 = s2*s2;
R2 = rb[2]+s*rb[3];s6 = s4*s2;
S2 = sb[2]+s*sb[3];
R3 = rb[4]+s*rb[5];
S3 = sb[4]+s*sb[5];
S4 = sb[6]+s*sb[7];
R = R1 + s2*R2 + s4*R3 + s6*rb[6];
S = S1 + s2*S2 + s4*S3 + s6*S4;
}
z = x;
SET_LOW_WORD(z,0);
r = __ieee754_exp(-z*z-0.5625)*
__ieee754_exp((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x;
} else {
if(hx>0) return tiny*tiny; else return two-tiny;
}
if (ix < 0x3ff40000) /* 0.84375 <= |x| < 1.25 */
{
double s2, s4, s6, P1, P2, P3, P4, Q1, Q2, Q3, Q4;
s = fabs (x) - one;
P1 = pa[0] + s * pa[1]; s2 = s * s;
Q1 = one + s * qa[1]; s4 = s2 * s2;
P2 = pa[2] + s * pa[3]; s6 = s4 * s2;
Q2 = qa[2] + s * qa[3];
P3 = pa[4] + s * pa[5];
Q3 = qa[4] + s * qa[5];
P4 = pa[6];
Q4 = qa[6];
P = P1 + s2 * P2 + s4 * P3 + s6 * P4;
Q = Q1 + s2 * Q2 + s4 * Q3 + s6 * Q4;
if (hx >= 0)
{
z = one - erx; return z - P / Q;
}
else
{
z = erx + P / Q; return one + z;
}
}
if (ix < 0x403c0000) /* |x|<28 */
{
x = fabs (x);
s = one / (x * x);
if (ix < 0x4006DB6D) /* |x| < 1/.35 ~ 2.857143*/
{
double R1, R2, R3, R4, S1, S2, S3, S4, s2, s4, s6, s8;
R1 = ra[0] + s * ra[1]; s2 = s * s;
S1 = one + s * sa[1]; s4 = s2 * s2;
R2 = ra[2] + s * ra[3]; s6 = s4 * s2;
S2 = sa[2] + s * sa[3]; s8 = s4 * s4;
R3 = ra[4] + s * ra[5];
S3 = sa[4] + s * sa[5];
R4 = ra[6] + s * ra[7];
S4 = sa[6] + s * sa[7];
R = R1 + s2 * R2 + s4 * R3 + s6 * R4;
S = S1 + s2 * S2 + s4 * S3 + s6 * S4 + s8 * sa[8];
}
else /* |x| >= 1/.35 ~ 2.857143 */
{
double R1, R2, R3, S1, S2, S3, S4, s2, s4, s6;
if (hx < 0 && ix >= 0x40180000)
return two - tiny; /* x < -6 */
R1 = rb[0] + s * rb[1]; s2 = s * s;
S1 = one + s * sb[1]; s4 = s2 * s2;
R2 = rb[2] + s * rb[3]; s6 = s4 * s2;
S2 = sb[2] + s * sb[3];
R3 = rb[4] + s * rb[5];
S3 = sb[4] + s * sb[5];
S4 = sb[6] + s * sb[7];
R = R1 + s2 * R2 + s4 * R3 + s6 * rb[6];
S = S1 + s2 * S2 + s4 * S3 + s6 * S4;
}
z = x;
SET_LOW_WORD (z, 0);
r = __ieee754_exp (-z * z - 0.5625) *
__ieee754_exp ((z - x) * (z + x) + R / S);
if (hx > 0)
return r / x;
else
return two - r / x;
}
else
{
if (hx > 0)
return tiny * tiny;
else
return two - tiny;
}
}
weak_alias (__erfc, erfc)
#ifdef NO_LONG_DOUBLE

223
sysdeps/ieee754/dbl-64/s_expm1.c
View File

@ -11,7 +11,7 @@
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
for performance improvement on pipelined processors.
*/
*/
/* expm1(x)
* Returns exp(x)-1, the exponential of x minus 1.
@ -113,116 +113,145 @@
#include <math_private.h>
#define one Q[0]
static const double
huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */
invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */
/* scaled coefficients related to expm1 */
Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */
1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
-7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
-2.01099218183624371326e-07}; /* BE8AFDB7 6E09C32D */
huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
ln2_hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
/* scaled coefficients related to expm1 */
Q[] = { 1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */
1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
-7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
-2.01099218183624371326e-07 }; /* BE8AFDB7 6E09C32D */
double
__expm1(double x)
__expm1 (double x)
{
double y,hi,lo,c,t,e,hxs,hfx,r1,h2,h4,R1,R2,R3;
int32_t k,xsb;
u_int32_t hx;
double y, hi, lo, c, t, e, hxs, hfx, r1, h2, h4, R1, R2, R3;
int32_t k, xsb;
u_int32_t hx;
GET_HIGH_WORD(hx,x);
xsb = hx&0x80000000; /* sign bit of x */
if(xsb==0) y=x; else y= -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */
GET_HIGH_WORD (hx, x);
xsb = hx & 0x80000000; /* sign bit of x */
if (xsb == 0)
y = x;
else
y = -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out huge and non-finite argument */
if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
u_int32_t low;
GET_LOW_WORD(low,x);
if(((hx&0xfffff)|low)!=0)
return x+x; /* NaN */
else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
}
if(x > o_threshold) {
__set_errno (ERANGE);
return huge*huge; /* overflow */
}
/* filter out huge and non-finite argument */
if (hx >= 0x4043687A) /* if |x|>=56*ln2 */
{
if (hx >= 0x40862E42) /* if |x|>=709.78... */
{
if (hx >= 0x7ff00000)
{
u_int32_t low;
GET_LOW_WORD (low, x);
if (((hx & 0xfffff) | low) != 0)
return x + x; /* NaN */
else
return (xsb == 0) ? x : -1.0; /* exp(+-inf)={inf,-1} */
}
if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */
math_force_eval(x+tiny); /* raise inexact */
return tiny-one; /* return -1 */
if (x > o_threshold)
{
__set_errno (ERANGE);
return huge * huge; /* overflow */
}
}
if (xsb != 0) /* x < -56*ln2, return -1.0 with inexact */
{
math_force_eval (x + tiny); /* raise inexact */
return tiny - one; /* return -1 */
}
}
/* argument reduction */
if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
if(xsb==0)
{hi = x - ln2_hi; lo = ln2_lo; k = 1;}
else
{hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
} else {
k = invln2*x+((xsb==0)?0.5:-0.5);
t = k;
hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
lo = t*ln2_lo;
/* argument reduction */
if (hx > 0x3fd62e42) /* if |x| > 0.5 ln2 */
{
if (hx < 0x3FF0A2B2) /* and |x| < 1.5 ln2 */
{
if (xsb == 0)
{
hi = x - ln2_hi; lo = ln2_lo; k = 1;
}
else
{
hi = x + ln2_hi; lo = -ln2_lo; k = -1;
}
x = hi - lo;
c = (hi-x)-lo;
}
else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */
t = huge+x; /* return x with inexact flags when x!=0 */
return x - (t-(huge+x));
else
{
k = invln2 * x + ((xsb == 0) ? 0.5 : -0.5);
t = k;
hi = x - t * ln2_hi; /* t*ln2_hi is exact here */
lo = t * ln2_lo;
}
else k = 0;
x = hi - lo;
c = (hi - x) - lo;
}
else if (hx < 0x3c900000) /* when |x|<2**-54, return x */
{
t = huge + x; /* return x with inexact flags when x!=0 */
return x - (t - (huge + x));
}
else
k = 0;
/* x is now in primary range */
hfx = 0.5*x;
hxs = x*hfx;
R1 = one+hxs*Q[1]; h2 = hxs*hxs;
R2 = Q[2]+hxs*Q[3]; h4 = h2*h2;
R3 = Q[4]+hxs*Q[5];
r1 = R1 + h2*R2 + h4*R3;
t = 3.0-r1*hfx;
e = hxs*((r1-t)/(6.0 - x*t));
if(k==0) return x - (x*e-hxs); /* c is 0 */
else {
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return 0.5*(x-e)-0.5;
if(k==1) {
if(x < -0.25) return -2.0*(e-(x+0.5));
else return one+2.0*(x-e);
}
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
u_int32_t high;
y = one-(e-x);
GET_HIGH_WORD(high,y);
SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */
return y-one;
}
t = one;
if(k<20) {
u_int32_t high;
SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
y = t-(e-x);
GET_HIGH_WORD(high,y);
SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */
} else {
u_int32_t high;
SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */
y = x-(e+t);
y += one;
GET_HIGH_WORD(high,y);
SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */
}
/* x is now in primary range */
hfx = 0.5 * x;
hxs = x * hfx;
R1 = one + hxs * Q[1]; h2 = hxs * hxs;
R2 = Q[2] + hxs * Q[3]; h4 = h2 * h2;
R3 = Q[4] + hxs * Q[5];
r1 = R1 + h2 * R2 + h4 * R3;
t = 3.0 - r1 * hfx;
e = hxs * ((r1 - t) / (6.0 - x * t));
if (k == 0)
return x - (x * e - hxs); /* c is 0 */
else
{
e = (x * (e - c) - c);
e -= hxs;
if (k == -1)
return 0.5 * (x - e) - 0.5;
if (k == 1)
{
if (x < -0.25)
return -2.0 * (e - (x + 0.5));
else
return one + 2.0 * (x - e);
}
return y;
if (k <= -2 || k > 56) /* suffice to return exp(x)-1 */
{
u_int32_t high;
y = one - (e - x);
GET_HIGH_WORD (high, y);
SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */
return y - one;
}
t = one;
if (k < 20)
{
u_int32_t high;
SET_HIGH_WORD (t, 0x3ff00000 - (0x200000 >> k)); /* t=1-2^-k */
y = t - (e - x);
GET_HIGH_WORD (high, y);
SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */
}
else
{
u_int32_t high;
SET_HIGH_WORD (t, ((0x3ff - k) << 20)); /* 2^-k */
y = x - (e + t);
y += one;
GET_HIGH_WORD (high, y);
SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */
}
}
return y;
}
weak_alias (__expm1, expm1)
#ifdef NO_LONG_DOUBLE

11
sysdeps/ieee754/dbl-64/s_fabs.c
View File

@ -21,12 +21,13 @@ static char rcsid[] = "$NetBSD: s_fabs.c,v 1.7 1995/05/10 20:47:13 jtc Exp $";
#include <math.h>
#include <math_private.h>
double __fabs(double x)
double
__fabs (double x)
{
u_int32_t high;
GET_HIGH_WORD(high,x);
SET_HIGH_WORD(x,high&0x7fffffff);
return x;
u_int32_t high;
GET_HIGH_WORD (high, x);
SET_HIGH_WORD (x, high & 0x7fffffff);
return x;
}
weak_alias (__fabs, fabs)
#ifdef NO_LONG_DOUBLE

9
sysdeps/ieee754/dbl-64/s_finite.c
View File

@ -23,11 +23,12 @@ static char rcsid[] = "$NetBSD: s_finite.c,v 1.8 1995/05/10 20:47:17 jtc Exp $";
#include <math_private.h>
#undef __finite
int __finite(double x)
int
__finite (double x)
{
int32_t hx;
GET_HIGH_WORD(hx,x);
return (int)((u_int32_t)((hx&0x7fffffff)-0x7ff00000)>>31);
int32_t hx;
GET_HIGH_WORD (hx, x);
return (int) ((u_int32_t) ((hx & 0x7fffffff) - 0x7ff00000) >> 31);
}
hidden_def (__finite)
weak_alias (__finite, finite)

91
sysdeps/ieee754/dbl-64/s_floor.c
View File

@ -24,44 +24,67 @@
static const double huge = 1.0e300;
double __floor(double x)
double
__floor (double x)
{
int32_t i0,i1,j0;
u_int32_t i,j;
EXTRACT_WORDS(i0,i1,x);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) { /* raise inexact if x != 0 */
math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */
if(i0>=0) {i0=i1=0;}
else if(((i0&0x7fffffff)|i1)!=0)
{ i0=0xbff00000;i1=0;}
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) return x; /* x is integral */
math_force_eval(huge+x); /* raise inexact flag */
if(i0<0) i0 += (0x00100000)>>j0;
i0 &= (~i); i1=0;
int32_t i0, i1, j0;
u_int32_t i, j;
EXTRACT_WORDS (i0, i1, x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20)
{
if (j0 < 0) /* raise inexact if x != 0 */
{
math_force_eval (huge + x); /* return 0*sign(x) if |x|<1 */
if (i0 >= 0)
{
i0 = i1 = 0;
}
} else if (j0>51) {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) return x; /* x is integral */
math_force_eval(huge+x); /* raise inexact flag */
if(i0<0) {
if(j0==20) i0+=1;
else {
j = i1+(1<<(52-j0));
if(j<i1) i0 +=1 ; /* got a carry */
i1=j;
}
else if (((i0 & 0x7fffffff) | i1) != 0)
{
i0 = 0xbff00000; i1 = 0;
}
i1 &= (~i);
}
INSERT_WORDS(x,i0,i1);
return x;
else
{
i = (0x000fffff) >> j0;
if (((i0 & i) | i1) == 0)
return x; /* x is integral */
math_force_eval (huge + x); /* raise inexact flag */
if (i0 < 0)
i0 += (0x00100000) >> j0;
i0 &= (~i); i1 = 0;
}
}
else if (j0 > 51)
{
if (j0 == 0x400)
return x + x; /* inf or NaN */
else
return x; /* x is integral */
}
else
{
i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0)
return x; /* x is integral */
math_force_eval (huge + x); /* raise inexact flag */
if (i0 < 0)
{
if (j0 == 20)
i0 += 1;
else
{
j = i1 + (1 << (52 - j0));
if (j < i1)
i0 += 1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
INSERT_WORDS (x, i0, i1);
return x;
}
#ifndef __floor
weak_alias (__floor, floor)

37
sysdeps/ieee754/dbl-64/s_frexp.c
View File

@ -28,25 +28,28 @@ static char rcsid[] = "$NetBSD: s_frexp.c,v 1.9 1995/05/10 20:47:24 jtc Exp $";
#include <math_private.h>
static const double
two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
double __frexp(double x, int *eptr)
double
__frexp (double x, int *eptr)
{
int32_t hx, ix, lx;
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
*eptr = 0;
if(ix>=0x7ff00000||((ix|lx)==0)) return x; /* 0,inf,nan */
if (ix<0x00100000) { /* subnormal */
x *= two54;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
*eptr = -54;
}
*eptr += (ix>>20)-1022;
hx = (hx&0x800fffff)|0x3fe00000;
SET_HIGH_WORD(x,hx);
return x;
int32_t hx, ix, lx;
EXTRACT_WORDS (hx, lx, x);
ix = 0x7fffffff & hx;
*eptr = 0;
if (ix >= 0x7ff00000 || ((ix | lx) == 0))
return x; /* 0,inf,nan */
if (ix < 0x00100000) /* subnormal */
{
x *= two54;
GET_HIGH_WORD (hx, x);
ix = hx & 0x7fffffff;
*eptr = -54;
}
*eptr += (ix >> 20) - 1022;
hx = (hx & 0x800fffff) | 0x3fe00000;
SET_HIGH_WORD (x, hx);
return x;
}
weak_alias (__frexp, frexp)
#ifdef NO_LONG_DOUBLE

10
sysdeps/ieee754/dbl-64/s_isinf.c
View File

@ -19,11 +19,11 @@ static char rcsid[] = "$NetBSD: s_isinf.c,v 1.3 1995/05/11 23:20:14 jtc Exp $";
int
__isinf (double x)
{
int32_t hx,lx;
EXTRACT_WORDS(hx,lx,x);
lx |= (hx & 0x7fffffff) ^ 0x7ff00000;
lx |= -lx;
return ~(lx >> 31) & (hx >> 30);
int32_t hx, lx;
EXTRACT_WORDS (hx, lx, x);
lx |= (hx & 0x7fffffff) ^ 0x7ff00000;
lx |= -lx;
return ~(lx >> 31) & (hx >> 30);
}
hidden_def (__isinf)
weak_alias (__isinf, isinf)

6
sysdeps/ieee754/dbl-64/s_isinf_ns.c
View File

@ -14,7 +14,7 @@
int
__isinf_ns (double x)
{
int32_t hx,lx;
EXTRACT_WORDS(hx,lx,x);
return !(lx | ((hx & 0x7fffffff) ^ 0x7ff00000));
int32_t hx, lx;
EXTRACT_WORDS (hx, lx, x);
return !(lx | ((hx & 0x7fffffff) ^ 0x7ff00000));
}

15
sysdeps/ieee754/dbl-64/s_isnan.c
View File

@ -23,14 +23,15 @@ static char rcsid[] = "$NetBSD: s_isnan.c,v 1.8 1995/05/10 20:47:36 jtc Exp $";
#include <math_private.h>
#undef __isnan
int __isnan(double x)
int
__isnan (double x)
{
int32_t hx,lx;
EXTRACT_WORDS(hx,lx,x);
hx &= 0x7fffffff;
hx |= (u_int32_t)(lx|(-lx))>>31;
hx = 0x7ff00000 - hx;
return (int)(((u_int32_t)hx)>>31);
int32_t hx, lx;
EXTRACT_WORDS (hx, lx, x);
hx &= 0x7fffffff;
hx |= (u_int32_t) (lx | (-lx)) >> 31;
hx = 0x7ff00000 - hx;
return (int) (((u_int32_t) hx) >> 31);
}
hidden_def (__isnan)
weak_alias (__isnan, isnan)

2
sysdeps/ieee754/dbl-64/s_llround.c
View File

@ -66,7 +66,7 @@ __llround (double x)
else
{
/* The number is too large. It is left implementation defined
what happens. */
what happens. */
return (long long int) x;
}

169
sysdeps/ieee754/dbl-64/s_log1p.c
View File

@ -11,7 +11,7 @@
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
for performance improvement on pipelined processors.
*/
*/
/* double log1p(double x)
*
@ -82,87 +82,112 @@
#include <math_private.h>
static const double
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lp[] = {0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */
3.999999999940941908e-01, /* 3FD99999 9997FA04 */
2.857142874366239149e-01, /* 3FD24924 94229359 */
2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
1.818357216161805012e-01, /* 3FC74664 96CB03DE */
1.531383769920937332e-01, /* 3FC39A09 D078C69F */
1.479819860511658591e-01}; /* 3FC2F112 DF3E5244 */
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lp[] = { 0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */
3.999999999940941908e-01, /* 3FD99999 9997FA04 */
2.857142874366239149e-01, /* 3FD24924 94229359 */
2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
1.818357216161805012e-01, /* 3FC74664 96CB03DE */
1.531383769920937332e-01, /* 3FC39A09 D078C69F */
1.479819860511658591e-01 }; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
double
__log1p(double x)
__log1p (double x)
{
double hfsq,f,c,s,z,R,u,z2,z4,z6,R1,R2,R3,R4;
int32_t k,hx,hu,ax;
double hfsq, f, c, s, z, R, u, z2, z4, z6, R1, R2, R3, R4;
int32_t k, hx, hu, ax;
GET_HIGH_WORD(hx,x);
ax = hx&0x7fffffff;
GET_HIGH_WORD (hx, x);
ax = hx & 0x7fffffff;
k = 1;
if (hx < 0x3FDA827A) { /* x < 0.41422 */
if(__builtin_expect(ax>=0x3ff00000, 0)) { /* x <= -1.0 */
if(x==-1.0) return -two54/(x-x);/* log1p(-1)=+inf */
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if(__builtin_expect(ax<0x3e200000, 0)) { /* |x| < 2**-29 */
math_force_eval(two54+x); /* raise inexact */
if (ax<0x3c900000) /* |x| < 2**-54 */
return x;
else
return x - x*x*0.5;
}
if(hx>0||hx<=((int32_t)0xbfd2bec3)) {
k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
k = 1;
if (hx < 0x3FDA827A) /* x < 0.41422 */
{
if (__builtin_expect (ax >= 0x3ff00000, 0)) /* x <= -1.0 */
{
if (x == -1.0)
return -two54 / (x - x); /* log1p(-1)=+inf */
else
return (x - x) / (x - x); /* log1p(x<-1)=NaN */
}
else if (__builtin_expect(hx >= 0x7ff00000, 0)) return x+x;
if(k!=0) {
if(hx<0x43400000) {
u = 1.0+x;
GET_HIGH_WORD(hu,u);
k = (hu>>20)-1023;
c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
c /= u;
} else {
u = x;
GET_HIGH_WORD(hu,u);
k = (hu>>20)-1023;
c = 0;
}
hu &= 0x000fffff;
if(hu<0x6a09e) {
SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */
} else {
k += 1;
SET_HIGH_WORD(u,hu|0x3fe00000); /* normalize u/2 */
hu = (0x00100000-hu)>>2;
}
f = u-1.0;
if (__builtin_expect (ax < 0x3e200000, 0)) /* |x| < 2**-29 */
{
math_force_eval (two54 + x); /* raise inexact */
if (ax < 0x3c900000) /* |x| < 2**-54 */
return x;
else
return x - x * x * 0.5;
}
hfsq=0.5*f*f;
if(hu==0) { /* |f| < 2**-20 */
if(f==zero) {
if(k==0) return zero;
else {c += k*ln2_lo; return k*ln2_hi+c;}
}
R = hfsq*(1.0-0.66666666666666666*f);
if(k==0) return f-R; else
return k*ln2_hi-((R-(k*ln2_lo+c))-f);
if (hx > 0 || hx <= ((int32_t) 0xbfd2bec3))
{
k = 0; f = x; hu = 1;
} /* -0.2929<x<0.41422 */
}
else if (__builtin_expect (hx >= 0x7ff00000, 0))
return x + x;
if (k != 0)
{
if (hx < 0x43400000)
{
u = 1.0 + x;
GET_HIGH_WORD (hu, u);
k = (hu >> 20) - 1023;
c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0); /* correction term */
c /= u;
}
s = f/(2.0+f);
z = s*s;
R1 = z*Lp[1]; z2=z*z;
R2 = Lp[2]+z*Lp[3]; z4=z2*z2;
R3 = Lp[4]+z*Lp[5]; z6=z4*z2;
R4 = Lp[6]+z*Lp[7];
R = R1 + z2*R2 + z4*R3 + z6*R4;
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
else
{
u = x;
GET_HIGH_WORD (hu, u);
k = (hu >> 20) - 1023;
c = 0;
}
hu &= 0x000fffff;
if (hu < 0x6a09e)
{
SET_HIGH_WORD (u, hu | 0x3ff00000); /* normalize u */
}
else
{
k += 1;
SET_HIGH_WORD (u, hu | 0x3fe00000); /* normalize u/2 */
hu = (0x00100000 - hu) >> 2;
}
f = u - 1.0;
}
hfsq = 0.5 * f * f;
if (hu == 0) /* |f| < 2**-20 */
{
if (f == zero)
{
if (k == 0)
return zero;
else
{
c += k * ln2_lo; return k * ln2_hi + c;
}
}
R = hfsq * (1.0 - 0.66666666666666666 * f);
if (k == 0)
return f - R;
else
return k * ln2_hi - ((R - (k * ln2_lo + c)) - f);
}
s = f / (2.0 + f);
z = s * s;
R1 = z * Lp[1]; z2 = z * z;
R2 = Lp[2] + z * Lp[3]; z4 = z2 * z2;
R3 = Lp[4] + z * Lp[5]; z6 = z4 * z2;
R4 = Lp[6] + z * Lp[7];
R = R1 + z2 * R2 + z4 * R3 + z6 * R4;
if (k == 0)
return f - (hfsq - s * (hfsq + R));
else
return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f);
}
weak_alias (__log1p, log1p)
#ifdef NO_LONG_DOUBLE

2
sysdeps/ieee754/dbl-64/s_logb.c
View File

@ -25,7 +25,7 @@ __logb (double x)
int32_t lx, ix, rix;
EXTRACT_WORDS (ix, lx, x);
ix &= 0x7fffffff; /* high |x| */
ix &= 0x7fffffff; /* high |x| */
if ((ix | lx) == 0)
return -1.0 / fabs (x);
if (ix >= 0x7ff00000)

2
sysdeps/ieee754/dbl-64/s_lrint.c
View File

@ -33,7 +33,7 @@ long int
__lrint (double x)
{
int32_t j0;
u_int32_t i0,i1;
u_int32_t i0, i1;
volatile double w;
double t;
long int result;

84
sysdeps/ieee754/dbl-64/s_modf.c
View File

@ -25,45 +25,59 @@
static const double one = 1.0;
double
__modf(double x, double *iptr)
__modf (double x, double *iptr)
{
int32_t i0,i1,j0;
u_int32_t i;
EXTRACT_WORDS(i0,i1,x);
j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */
if(j0<20) { /* integer part in high x */
if(j0<0) { /* |x|<1 */
INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */
return x;
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) { /* x is integral */
*iptr = x;
INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
return x;
} else {
INSERT_WORDS(*iptr,i0&(~i),0);
return x - *iptr;
}
int32_t i0, i1, j0;
u_int32_t i;
EXTRACT_WORDS (i0, i1, x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; /* exponent of x */
if (j0 < 20) /* integer part in high x */
{
if (j0 < 0) /* |x|<1 */
{
INSERT_WORDS (*iptr, i0 & 0x80000000, 0); /* *iptr = +-0 */
return x;
}
else
{
i = (0x000fffff) >> j0;
if (((i0 & i) | i1) == 0) /* x is integral */
{
*iptr = x;
INSERT_WORDS (x, i0 & 0x80000000, 0); /* return +-0 */
return x;
}
} else if (__builtin_expect(j0>51, 0)) { /* no fraction part */
*iptr = x*one;
/* We must handle NaNs separately. */
if (j0 == 0x400 && ((i0 & 0xfffff) | i1))
return x*one;
INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
return x;
} else { /* fraction part in low x */
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) { /* x is integral */
*iptr = x;
INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
return x;
} else {
INSERT_WORDS(*iptr,i0,i1&(~i));
return x - *iptr;
else
{
INSERT_WORDS (*iptr, i0 & (~i), 0);
return x - *iptr;
}
}
}
else if (__builtin_expect (j0 > 51, 0)) /* no fraction part */
{
*iptr = x * one;
/* We must handle NaNs separately. */
if (j0 == 0x400 && ((i0 & 0xfffff) | i1))
return x * one;
INSERT_WORDS (x, i0 & 0x80000000, 0); /* return +-0 */
return x;
}
else /* fraction part in low x */
{
i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0) /* x is integral */
{
*iptr = x;
INSERT_WORDS (x, i0 & 0x80000000, 0); /* return +-0 */
return x;
}
else
{
INSERT_WORDS (*iptr, i0, i1 & (~i));
return x - *iptr;
}
}
}
weak_alias (__modf, modf)
#ifdef NO_LONG_DOUBLE

63
sysdeps/ieee754/dbl-64/s_nearbyint.c
View File

@ -29,40 +29,47 @@ static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
#include <math_private.h>
static const double
TWO52[2]={
TWO52[2] = {
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
};
double __nearbyint(double x)
double
__nearbyint (double x)
{
fenv_t env;
int32_t i0,j0,sx;
double w,t;
GET_HIGH_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<52) {
if(j0<0) {
libc_feholdexcept (&env);
w = TWO52[sx]+x;
t = w-TWO52[sx];
math_force_eval (t);
libc_fesetenv (&env);
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
}
} else {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
fenv_t env;
int32_t i0, j0, sx;
double w, t;
GET_HIGH_WORD (i0, x);
sx = (i0 >> 31) & 1;
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 52)
{
if (j0 < 0)
{
libc_feholdexcept (&env);
w = TWO52[sx] + x;
t = w - TWO52[sx];
math_force_eval (t);
libc_fesetenv (&env);
GET_HIGH_WORD (i0, t);
SET_HIGH_WORD (t, (i0 & 0x7fffffff) | (sx << 31));
return t;
}
libc_feholdexcept (&env);
w = TWO52[sx]+x;
t = w-TWO52[sx];
math_force_eval (t);
libc_fesetenv (&env);
return t;
}
else
{
if (j0 == 0x400)
return x + x; /* inf or NaN */
else
return x; /* x is integral */
}
libc_feholdexcept (&env);
w = TWO52[sx] + x;
t = w - TWO52[sx];
math_force_eval (t);
libc_fesetenv (&env);
return t;
}
weak_alias (__nearbyint, nearbyint)
#ifdef NO_LONG_DOUBLE

16
sysdeps/ieee754/dbl-64/s_remquo.c
View File

@ -28,8 +28,8 @@ static const double zero = 0.0;
double
__remquo (double x, double y, int *quo)
{
int32_t hx,hy;
u_int32_t sx,lx,ly;
int32_t hx, hy;
u_int32_t sx, lx, ly;
int cquo, qs;
EXTRACT_WORDS (hx, lx, x);
@ -41,14 +41,14 @@ __remquo (double x, double y, int *quo)
/* Purge off exception values. */
if ((hy | ly) == 0)
return (x * y) / (x * y); /* y = 0 */
if ((hx >= 0x7ff00000) /* x not finite */
|| ((hy >= 0x7ff00000) /* p is NaN */
return (x * y) / (x * y); /* y = 0 */
if ((hx >= 0x7ff00000) /* x not finite */
|| ((hy >= 0x7ff00000) /* p is NaN */
&& (((hy - 0x7ff00000) | ly) != 0)))
return (x * y) / (x * y);
if (hy <= 0x7fbfffff)
x = __ieee754_fmod (x, 8 * y); /* now x < 8y */
x = __ieee754_fmod (x, 8 * y); /* now x < 8y */
if (((hx - hy) | (lx - ly)) == 0)
{
@ -56,8 +56,8 @@ __remquo (double x, double y, int *quo)
return zero * x;
}
x = fabs (x);
y = fabs (y);
x = fabs (x);
y = fabs (y);
cquo = 0;
if (x >= 4 * y)

46
sysdeps/ieee754/dbl-64/s_rint.c
View File

@ -24,33 +24,39 @@
#include <math_private.h>
static const double
TWO52[2]={
TWO52[2] = {
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
};
double
__rint(double x)
__rint (double x)
{
int32_t i0,j0,sx;
double w,t;
GET_HIGH_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<52) {
if(j0<0) {
w = TWO52[sx]+x;
t = w-TWO52[sx];
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
}
} else {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
int32_t i0, j0, sx;
double w, t;
GET_HIGH_WORD (i0, x);
sx = (i0 >> 31) & 1;
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 52)
{
if (j0 < 0)
{
w = TWO52[sx] + x;
t = w - TWO52[sx];
GET_HIGH_WORD (i0, t);
SET_HIGH_WORD (t, (i0 & 0x7fffffff) | (sx << 31));
return t;
}
w = TWO52[sx]+x;
return w-TWO52[sx];
}
else
{
if (j0 == 0x400)
return x + x; /* inf or NaN */
else
return x; /* x is integral */
}
w = TWO52[sx] + x;
return w - TWO52[sx];
}
#ifndef __rint
weak_alias (__rint, rint)

61
sysdeps/ieee754/dbl-64/s_scalbln.c
View File

@ -20,38 +20,43 @@
#include <math_private.h>
static const double
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
double
__scalbln (double x, long int n)
{
int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (__builtin_expect(k==0, 0)) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx,x);
k = ((hx&0x7ff00000)>>20) - 54;
}
if (__builtin_expect(k==0x7ff, 0)) return x+x; /* NaN or Inf */
if (__builtin_expect(n< -50000, 0))
return tiny*__copysign(tiny,x); /*underflow*/
if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0))
return huge*__copysign(huge,x); /* overflow */
/* Now k and n are bounded we know that k = k+n does not
overflow. */
k = k+n;
if (__builtin_expect(k > 0, 1)) /* normal result */
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
if (k <= -54)
return tiny*__copysign(tiny,x); /*underflow*/
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;
int32_t k, hx, lx;
EXTRACT_WORDS (hx, lx, x);
k = (hx & 0x7ff00000) >> 20; /* extract exponent */
if (__builtin_expect (k == 0, 0)) /* 0 or subnormal x */
{
if ((lx | (hx & 0x7fffffff)) == 0)
return x; /* +-0 */
x *= two54;
GET_HIGH_WORD (hx, x);
k = ((hx & 0x7ff00000) >> 20) - 54;
}
if (__builtin_expect (k == 0x7ff, 0))
return x + x; /* NaN or Inf */
if (__builtin_expect (n < -50000, 0))
return tiny * __copysign (tiny, x); /*underflow*/
if (__builtin_expect (n > 50000 || k + n > 0x7fe, 0))
return huge * __copysign (huge, x); /* overflow */
/* Now k and n are bounded we know that k = k+n does not
overflow. */
k = k + n;
if (__builtin_expect (k > 0, 1)) /* normal result */
{
SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x;
}
if (k <= -54)
return tiny * __copysign (tiny, x); /*underflow*/
k += 54; /* subnormal result */
SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20));
return x * twom54;
}
weak_alias (__scalbln, scalbln)
#ifdef NO_LONG_DOUBLE

61
sysdeps/ieee754/dbl-64/s_scalbn.c
View File

@ -20,38 +20,43 @@
#include <math_private.h>
static const double
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
double
__scalbn (double x, int n)
{
int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (__builtin_expect(k==0, 0)) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx,x);
k = ((hx&0x7ff00000)>>20) - 54;
}
if (__builtin_expect(k==0x7ff, 0)) return x+x; /* NaN or Inf */
if (__builtin_expect(n< -50000, 0))
return tiny*__copysign(tiny,x); /*underflow*/
if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0))
return huge*__copysign(huge,x); /* overflow */
/* Now k and n are bounded we know that k = k+n does not
overflow. */
k = k+n;
if (__builtin_expect(k > 0, 1)) /* normal result */
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
if (k <= -54)
return tiny*__copysign(tiny,x); /*underflow*/
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;
int32_t k, hx, lx;
EXTRACT_WORDS (hx, lx, x);
k = (hx & 0x7ff00000) >> 20; /* extract exponent */
if (__builtin_expect (k == 0, 0)) /* 0 or subnormal x */
{
if ((lx | (hx & 0x7fffffff)) == 0)
return x; /* +-0 */
x *= two54;
GET_HIGH_WORD (hx, x);
k = ((hx & 0x7ff00000) >> 20) - 54;
}
if (__builtin_expect (k == 0x7ff, 0))
return x + x; /* NaN or Inf */
if (__builtin_expect (n < -50000, 0))
return tiny * __copysign (tiny, x); /*underflow*/
if (__builtin_expect (n > 50000 || k + n > 0x7fe, 0))
return huge * __copysign (huge, x); /* overflow */
/* Now k and n are bounded we know that k = k+n does not
overflow. */
k = k + n;
if (__builtin_expect (k > 0, 1)) /* normal result */
{
SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20)); return x;
}
if (k <= -54)
return tiny * __copysign (tiny, x); /*underflow*/
k += 54; /* subnormal result */
SET_HIGH_WORD (x, (hx & 0x800fffff) | (k << 20));
return x * twom54;
}
weak_alias (__scalbn, scalbn)
#ifdef NO_LONG_DOUBLE

2
sysdeps/ieee754/dbl-64/s_sin.c
View File

@ -725,7 +725,7 @@ slow (double x)
}
/*******************************************************************************/
/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
/* and if result still doesn't accurate enough by mpsin or dubsin */
/*******************************************************************************/

2
sysdeps/ieee754/dbl-64/s_sincos.c
View File

@ -32,7 +32,7 @@ __sincos (double x, double *sinx, double *cosx)
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if (ix>=0x7ff00000)
if (ix >= 0x7ff00000)
{
/* sin(Inf or NaN) is NaN */
*sinx = *cosx = x - x;

8
sysdeps/ieee754/dbl-64/s_tan.c
View File

@ -59,8 +59,8 @@ tan (double x)
int ux, i, n;
double a, da, a2, b, db, c, dc, c1, cc1, c2, cc2, c3, cc3, fi, ffi, gi, pz,
s, sy, t, t1, t2, t3, t4, t7, t8, t9, t10, w, x2, xn, xx2, y, ya, yya, z0,
z, zz, z2, zz2;
s, sy, t, t1, t2, t3, t4, t7, t8, t9, t10, w, x2, xn, xx2, y, ya,
yya, z0, z, zz, z2, zz2;
#ifndef DLA_FMS
double t5, t6;
#endif
@ -98,7 +98,6 @@ tan (double x)
/* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
if (w <= g2.d)
{
/* First stage */
x2 = x * x;
@ -150,7 +149,6 @@ tan (double x)
/* (III) The case 0.0608 < abs(x) <= 0.787 */
if (w <= g3.d)
{
/* First stage */
i = ((int) (mfftnhf.d + TWO8 * w));
z = w - xfg[i][0].d;
@ -377,7 +375,7 @@ tan (double x)
/* Second stage */
ffi = xfg[i][3].d;
EADD (z0, yya, z, zz)
MUL2 (z, zz, z, zz, z2, zz2, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (z, zz, z, zz, z2, zz2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);

68
sysdeps/ieee754/dbl-64/s_tanh.c
View File

@ -41,41 +41,51 @@ static char rcsid[] = "$NetBSD: s_tanh.c,v 1.7 1995/05/10 20:48:22 jtc Exp $";
#include <math.h>
#include <math_private.h>
static const double one=1.0, two=2.0, tiny = 1.0e-300;
static const double one = 1.0, two = 2.0, tiny = 1.0e-300;
double __tanh(double x)
double
__tanh (double x)
{
double t,z;
int32_t jx,ix,lx;
double t, z;
int32_t jx, ix, lx;
/* High word of |x|. */
EXTRACT_WORDS(jx,lx,x);
ix = jx&0x7fffffff;
/* High word of |x|. */
EXTRACT_WORDS (jx, lx, x);
ix = jx & 0x7fffffff;
/* x is INF or NaN */
if(ix>=0x7ff00000) {
if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */
else return one/x-one; /* tanh(NaN) = NaN */
/* x is INF or NaN */
if (ix >= 0x7ff00000)
{
if (jx >= 0)
return one / x + one; /* tanh(+-inf)=+-1 */
else
return one / x - one; /* tanh(NaN) = NaN */
}
/* |x| < 22 */
if (ix < 0x40360000) /* |x|<22 */
{
if ((ix | lx) == 0)
return x; /* x == +-0 */
if (ix < 0x3c800000) /* |x|<2**-55 */
return x * (one + x); /* tanh(small) = small */
if (ix >= 0x3ff00000) /* |x|>=1 */
{
t = __expm1 (two * fabs (x));
z = one - two / (t + two);
}
/* |x| < 22 */
if (ix < 0x40360000) { /* |x|<22 */
if ((ix | lx) == 0)
return x; /* x == +-0 */
if (ix<0x3c800000) /* |x|<2**-55 */
return x*(one+x); /* tanh(small) = small */
if (ix>=0x3ff00000) { /* |x|>=1 */
t = __expm1(two*fabs(x));
z = one - two/(t+two);
} else {
t = __expm1(-two*fabs(x));
z= -t/(t+two);
}
/* |x| > 22, return +-1 */
} else {
z = one - tiny; /* raised inexact flag */
else
{
t = __expm1 (-two * fabs (x));
z = -t / (t + two);
}
return (jx>=0)? z: -z;
/* |x| > 22, return +-1 */
}
else
{
z = one - tiny; /* raised inexact flag */
}
return (jx >= 0) ? z : -z;
}
weak_alias (__tanh, tanh)
#ifdef NO_LONG_DOUBLE