gered/picoc
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity
picoc/math_library.c
zik.saleeba a38afda0ab Limited decimal printing to sane precision. 
Made C library able to access new math library.
UNIX compile defaults to using new math library.
Changed name of new math library to math.library.c
Added defines for new math library.
Renamed new soft floating point library to softfloat_library.c

git-svn-id: http://picoc.googlecode.com/svn/trunk@304 21eae674-98b7-11dd-bd71-f92a316d2d60
2009-05-29 00:12:12 +00:00

2829 lines
80 KiB
C
Raw Blame History

/*
* This is a modified version of FPLIBM for picoc.
*
* FDLIBM (Freely Distributable LIBM) is a C math library
* for machines that support IEEE 754 floating-point arithmetic.
* In this release, only double precision is supported.
*
* FDLIBM is intended to provide a reasonably portable (see
* assumptions below), reference quality (below one ulp for
* major functions like sin,cos,exp,log) math library
* (libm.a). For a copy of FDLIBM, please see
* http://www.netlib.org/fdlibm/
* or
* http://www.validlab.com/software/
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include "picoc.h"
#ifdef NEED_MATH_LIBRARY
/* Sometimes it's necessary to define BIG_ENDIAN explicitly */
#ifndef BIG_ENDIAN
#define __HI(x) *(1+(int*)&x)
#define __LO(x) *(int*)&x
#define __HIp(x) *(1+(int*)x)
#define __LOp(x) *(int*)x
#else
#define __HI(x) *(int*)&x
#define __LO(x) *(1+(int*)&x)
#define __HIp(x) *(int*)x
#define __LOp(x) *(1+(int*)x)
#endif
/* handy constants */
static const double
tiny= 1.00000000000000000000e-300,
huge= 1.00000000000000000000e+300,
shuge = 1.0e307,
halF[2] = {0.5,-0.5,},
zero= 0.0,
one= 1.0, /* 0x3FF00000, 0x00000000 */
two = 2.0,
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08, /* 0x3E700000, 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 */
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
pio4 = 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
pio4lo = 3.06161699786838301793e-17, /* 3C81A626, 33145C07 */
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
qS4 = 7.70381505559019352791e-02, /* 0x3FB3B8C5, 0xB12E9282 */
twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
-6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
-1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01, /* 3FC2F112 DF3E5244 */
ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13, /* 0x3D59FEF3, 0x11F12B36 */
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
ivln2_l = 1.92596299112661746887e-08, /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10, /* 0x3DE5D93A, 0x5ACFD57C */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C6 = -1.13596475577881948265e-11, /* 0xBDA8FAE9, 0xBE8838D4 */
T[] = {
3.33333333333334091986e-01, /* 3FD55555, 55555563 */
1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */
5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */
2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */
8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */
3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */
1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */
5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */
2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */
7.81794442939557092300e-05, /* 3F147E88, A03792A6 */
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
};
/*
* __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
* double x[],y[]; int e0,nx,prec; int ipio2[];
*
* __kernel_rem_pio2 return the last three digits of N with
* y = x - N*pi/2
* so that |y| < pi/2.
*
* The method is to compute the integer (mod 8) and fraction parts of
* (2/pi)*x without doing the full multiplication. In general we
* skip the part of the product that are known to be a huge integer (
* more accurately, = 0 mod 8 ). Thus the number of operations are
* independent of the exponent of the input.
*
* (2/pi) is represented by an array of 24-bit integers in ipio2[].
*
* Input parameters:
* x[] The input value (must be positive) is broken into nx
* pieces of 24-bit integers in double precision format.
* x[i] will be the i-th 24 bit of x. The scaled exponent
* of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
* match x's up to 24 bits.
*
* Example of breaking a double positive z into x[0]+x[1]+x[2]:
* e0 = ilogb(z)-23
* z = scalbn(z,-e0)
* for i = 0,1,2
* x[i] = floor(z)
* z = (z-x[i])*2**24
*
*
* y[] ouput result in an array of double precision numbers.
* The dimension of y[] is:
* 24-bit precision 1
* 53-bit precision 2
* 64-bit precision 2
* 113-bit precision 3
* The actual value is the sum of them. Thus for 113-bit
* precison, one may have to do something like:
*
* long double t,w,r_head, r_tail;
* t = (long double)y[2] + (long double)y[1];
* w = (long double)y[0];
* r_head = t+w;
* r_tail = w - (r_head - t);
*
* e0 The exponent of x[0]
*
* nx dimension of x[]
*
* prec an integer indicating the precision:
* 0 24 bits (single)
* 1 53 bits (double)
* 2 64 bits (extended)
* 3 113 bits (quad)
*
* ipio2[]
* integer array, contains the (24*i)-th to (24*i+23)-th
* bit of 2/pi after binary point. The corresponding
* floating value is
*
* ipio2[i] * 2^(-24(i+1)).
*
* External function:
* double scalbn(), floor();
*
*
* Here is the description of some local variables:
*
* jk jk+1 is the initial number of terms of ipio2[] needed
* in the computation. The recommended value is 2,3,4,
* 6 for single, double, extended,and quad.
*
* jz local integer variable indicating the number of
* terms of ipio2[] used.
*
* jx nx - 1
*
* jv index for pointing to the suitable ipio2[] for the
* computation. In general, we want
* ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
* is an integer. Thus
* e0-3-24*jv >= 0 or (e0-3)/24 >= jv
* Hence jv = max(0,(e0-3)/24).
*
* jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
*
* q[] double array with integral value, representing the
* 24-bits chunk of the product of x and 2/pi.
*
* q0 the corresponding exponent of q[0]. Note that the
* exponent for q[i] would be q0-24*i.
*
* PIo2[] double precision array, obtained by cutting pi/2
* into 24 bits chunks.
*
* f[] ipio2[] in floating point
*
* iq[] integer array by breaking up q[] in 24-bits chunk.
*
* fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
*
* ih integer. If >0 it indicates q[] is >= 0.5, hence
* it also indicates the *sign* of the result.
*
*/
/*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
static const double PIo2[] = {
1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
};
int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2)
{
int 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;
/* 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];
/* 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;
recompute:
/* distill q[] into iq[] reversingly */
for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
fw = (double)((int)(twon24* z));
iq[i] = (int)(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 = (int) 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;
}
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)((int)(twon24*z));
iq[jz] = (int)(z-two24*fw);
jz += 1; q0 += 24;
iq[jz] = (int) fw;
} else iq[jz] = (int) 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:
fw = 0.0;
for (i=jz;i>=0;i--) fw += fq[i];
y[0] = (ih==0)? fw: -fw;
fw = fq[0]-fw;
for (i=1;i<=jz;i++) fw += fq[i];
y[1] = (ih==0)? fw: -fw;
break;
case 3: /* painful */
for (i=jz;i>0;i--) {
fw = fq[i-1]+fq[i];
fq[i] += fq[i-1]-fw;
fq[i-1] = fw;
}
for (i=jz;i>1;i--) {
fw = fq[i-1]+fq[i];
fq[i] += fq[i-1]-fw;
fq[i-1] = fw;
}
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;
}
/* __ieee754_rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __kernel_rem_pio2()
*/
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
static const int two_over_pi[] = {
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
static const int npio2_hw[] = {
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
0x404858EB, 0x404921FB,
};
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 33 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 33 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
int __ieee754_rem_pio2(double x, double *y)
{
double z,w,t,r,fn;
double tx[3];
int e0,i,j,nx,n,ix,hx;
hx = __HI(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 = (int) (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 {
j = ix>>20;
y[0] = r-w;
i = j-(((__HI(y[0]))>>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;
i = j-(((__HI(y[0]))>>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;
}
/*
* 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) */
__LO(z) = __LO(x);
e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
__HI(z) = ix - (e0<<20);
for(i=0;i<2;i++) {
tx[i] = (double)((int)(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;
}
/* __ieee754_exp(x)
* Returns the exponential of x.
*
* Method
* 1. Argument reduction:
* Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
* Given x, find r and integer k such that
*
* x = k*ln2 + r, |r| <= 0.5*ln2.
*
* Here r will be represented as r = hi-lo for better
* accuracy.
*
* 2. Approximation of exp(r) by a special rational function on
* the interval [0,0.34658]:
* Write
* R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
* We use a special Remes algorithm on [0,0.34658] to generate
* a polynomial of degree 5 to approximate R. The maximum error
* of this polynomial approximation is bounded by 2**-59. In
* other words,
* R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
* (where z=r*r, and the values of P1 to P5 are listed below)
* and
* | 5 | -59
* | 2.0+P1*z+...+P5*z - R(z) | <= 2
* | |
* The computation of exp(r) thus becomes
* 2*r
* exp(r) = 1 + -------
* R - r
* r*R1(r)
* = 1 + r + ----------- (for better accuracy)
* 2 - R1(r)
* where
* 2 4 10
* R1(r) = r - (P1*r + P2*r + ... + P5*r ).
*
* 3. Scale back to obtain exp(x):
* From step 1, we have
* exp(x) = 2^k * exp(r)
*
* Special cases:
* exp(INF) is INF, exp(NaN) is NaN;
* exp(-INF) is 0, and
* for finite argument, only exp(0)=1 is exact.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Misc. info.
* For IEEE double
* if x > 7.09782712893383973096e+02 then exp(x) overflow
* if x < -7.45133219101941108420e+02 then exp(x) underflow
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
double __ieee754_exp(double x) /* default IEEE double exp */
{
double y,hi,lo,c,t;
int k,xsb;
unsigned hx;
hx = __HI(x); /* high word of x */
xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
if(((hx&0xfffff)|__LO(x))!=0)
return x+x; /* NaN */
else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
}
if(x > o_threshold) return huge*huge; /* overflow */
if(x < u_threshold) return twom1000*twom1000; /* underflow */
}
/* argument reduction */
if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
} else {
k = (int)(invln2*x+halF[xsb]);
t = k;
hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
lo = t*ln2LO[0];
}
x = hi - lo;
}
else if(hx < 0x3e300000) { /* when |x|<2**-28 */
if(huge+x>one) return one+x;/* trigger inexact */
}
else k = 0;
/* x is now in primary range */
t = x*x;
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);
if(k >= -1021) {
__HI(y) += (k<<20); /* add k to y's exponent */
return y;
} else {
__HI(y) += ((k+1000)<<20);/* add k to y's exponent */
return y*twom1000;
}
}
/* __ieee754_log(x)
* Return the logrithm of x
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
double __ieee754_log(double x)
{
double hfsq,f,s,z,R,w,t1,t2,dk;
int k,hx,i,j;
unsigned lx;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
__HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
k += (i>>20);
f = x-1.0;
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
if(f==zero) { if(k==0) return zero; else {dk=(double)k;
return dk*ln2_hi+dk*ln2_lo;} }
R = f*f*(0.5-0.33333333333333333*f);
if(k==0) return f-R; else {dk=(double)k;
return dk*ln2_hi-((R-dk*ln2_lo)-f);}
}
s = f/(2.0+f);
dk = (double)k;
z = s*s;
i = hx-0x6147a;
w = z*z;
j = 0x6b851-hx;
t1= w*(Lg2+w*(Lg4+w*Lg6));
t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
i |= j;
R = t2+t1;
if(i>0) {
hfsq=0.5*f*f;
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
} else {
if(k==0) return f-s*(f-R); else
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}
/* __kernel_sin( x, y, iy)
* kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
*
* Algorithm
* 1. Since sin(-x) = -sin(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
* 3. sin(x) is approximated by a polynomial of degree 13 on
* [0,pi/4]
* 3 13
* sin(x) ~ x + S1*x + ... + S6*x
* where
*
* |sin(x) 2 4 6 8 10 12 | -58
* |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
* | x |
*
* 4. sin(x+y) = sin(x) + sin'(x')*y
* ~ sin(x) + (1-x*x/2)*y
* For better accuracy, let
* 3 2 2 2 2
* r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
* then 3 2
* sin(x) = x + (S1*x + (x *(r-y/2)+y))
*/
double __kernel_sin(double x, double y, int iy)
{
double z,r,v;
int ix;
ix = __HI(x)&0x7fffffff; /* high word of x */
if(ix<0x3e400000) /* |x| < 2**-27 */
{if((int)x==0) return x;} /* generate inexact */
z = x*x;
v = z*x;
r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
if(iy==0) return x+v*(S1+z*r);
else return x-((z*(half*y-v*r)-y)-v*S1);
}
/*
* __kernel_cos( x, y )
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
*
* Algorithm
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
* 3. cos(x) is approximated by a polynomial of degree 14 on
* [0,pi/4]
* 4 14
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
* where the remez error is
*
* | 2 4 6 8 10 12 14 | -58
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
* | |
*
* 4 6 8 10 12 14
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
* cos(x) = 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* ~ cos(x) - x*y,
* a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
* For better accuracy when x > 0.3, let qx = |x|/4 with
* the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
* Then
* cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
* Note that 1-qx and (x*x/2-qx) is EXACT here, and the
* magnitude of the latter is at least a quarter of x*x/2,
* thus, reducing the rounding error in the subtraction.
*/
double __kernel_cos(double x, double y)
{
double a,hz,z,r,qx;
int ix;
ix = __HI(x)&0x7fffffff; /* ix = |x|'s high word*/
if(ix<0x3e400000) { /* if x < 2**27 */
if(((int)x)==0) return one; /* generate inexact */
}
z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
if(ix < 0x3FD33333) /* if |x| < 0.3 */
return one - (0.5*z - (z*r - x*y));
else {
if(ix > 0x3fe90000) { /* x > 0.78125 */
qx = 0.28125;
} else {
__HI(qx) = ix-0x00200000; /* x/4 */
__LO(qx) = 0;
}
hz = 0.5*z-qx;
a = one-qx;
return a - (hz - (z*r-x*y));
}
}
/* __kernel_tan( x, y, k )
* kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
*
* Algorithm
* 1. Since tan(-x) = -tan(x), we need only to consider positive x.
* 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
* 3. tan(x) is approximated by a odd polynomial of degree 27 on
* [0,0.67434]
* 3 27
* tan(x) ~ x + T1*x + ... + T13*x
* where
*
* |tan(x) 2 4 26 | -59.2
* |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
* | x |
*
* Note: tan(x+y) = tan(x) + tan'(x)*y
* ~ tan(x) + (1+x*x)*y
* Therefore, for better accuracy in computing tan(x+y), let
* 3 2 2 2 2
* r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
* then
* 3 2
* tan(x+y) = x + (T1*x + (x *(r+y)+y))
*
* 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
* tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
* = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
*/
double __kernel_tan(double x, double y, int iy)
{
double z, r, v, w, s;
int ix, hx;
hx = __HI(x); /* high word of x */
ix = hx & 0x7fffffff; /* high word of |x| */
if (ix < 0x3e300000) { /* x < 2**-28 */
if ((int) x == 0) { /* generate inexact */
if (((ix | __LO(x)) | (iy + 1)) == 0)
return one / fabs(x);
else {
if (iy == 1)
return x;
else { /* compute -1 / (x+y) carefully */
double a, t;
z = w = x + y;
__LO(z) = 0;
v = y - (z - x);
t = a = -one / w;
__LO(t) = 0;
s = one + t * z;
return t + a * (s + t * v);
}
}
}
}
if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */
if (hx < 0) {
x = -x;
y = -y;
}
z = pio4 - x;
w = pio4lo - y;
x = z + w;
y = 0.0;
}
z = x * x;
w = z * z;
/*
* Break x^5*(T[1]+x^2*T[2]+...) into
* x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
* x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
*/
r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] +
w * T[11]))));
v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] +
w * T[12])))));
s = z * x;
r = y + z * (s * (r + v) + y);
r += T[0] * s;
w = x + r;
if (ix >= 0x3FE59428) {
v = (double) iy;
return (double) (1 - ((hx >> 30) & 2)) *
(v - 2.0 * (x - (w * w / (w + v) - r)));
}
if (iy == 1)
return w;
else {
/*
* if allow error up to 2 ulp, simply return
* -1.0 / (x+r) here
*/
/* compute -1.0 / (x+r) accurately */
double a, t;
z = w;
__LO(z) = 0;
v = r - (z - x); /* z+v = r+x */
t = a = -1.0 / w; /* a = -1.0/w */
__LO(t) = 0;
s = 1.0 + t * z;
return t + a * (s + t * v);
}
}
/* sin(x)
* Return sine function of x.
*
* kernel function:
* __kernel_sin ... sine function on [-pi/4,pi/4]
* __kernel_cos ... cose function on [-pi/4,pi/4]
* __ieee754_rem_pio2 ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
* in [-pi/4 , +pi/4], and let n = k mod 4.
* We have
*
* n sin(x) cos(x) tan(x)
* ----------------------------------------------------------
* 0 S C T
* 1 C -S -1/T
* 2 -S -C T
* 3 -C S -1/T
* ----------------------------------------------------------
*
* Special cases:
* Let trig be any of sin, cos, or tan.
* trig(+-INF) is NaN, with signals;
* trig(NaN) is that NaN;
*
* Accuracy:
* TRIG(x) returns trig(x) nearly rounded
*/
double math_sin(double x)
{
double y[2],z=0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0);
/* sin(Inf or NaN) is NaN */
else if (ix>=0x7ff00000) return x-x;
/* argument reduction needed */
else {
n = __ieee754_rem_pio2(x,y);
switch(n&3) {
case 0: return __kernel_sin(y[0],y[1],1);
case 1: return __kernel_cos(y[0],y[1]);
case 2: return -__kernel_sin(y[0],y[1],1);
default:
return -__kernel_cos(y[0],y[1]);
}
}
}
/* cos(x)
* Return cosine function of x.
*
* kernel function:
* __kernel_sin ... sine function on [-pi/4,pi/4]
* __kernel_cos ... cosine function on [-pi/4,pi/4]
* __ieee754_rem_pio2 ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
* in [-pi/4 , +pi/4], and let n = k mod 4.
* We have
*
* n sin(x) cos(x) tan(x)
* ----------------------------------------------------------
* 0 S C T
* 1 C -S -1/T
* 2 -S -C T
* 3 -C S -1/T
* ----------------------------------------------------------
*
* Special cases:
* Let trig be any of sin, cos, or tan.
* trig(+-INF) is NaN, with signals;
* trig(NaN) is that NaN;
*
* Accuracy:
* TRIG(x) returns trig(x) nearly rounded
*/
double math_cos(double x)
{
double y[2],z=0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
/* cos(Inf or NaN) is NaN */
else if (ix>=0x7ff00000) return x-x;
/* argument reduction needed */
else {
n = __ieee754_rem_pio2(x,y);
switch(n&3) {
case 0: return __kernel_cos(y[0],y[1]);
case 1: return -__kernel_sin(y[0],y[1],1);
case 2: return -__kernel_cos(y[0],y[1]);
default:
return __kernel_sin(y[0],y[1],1);
}
}
}
/* tan(x)
* Return tangent function of x.
*
* kernel function:
* __kernel_tan ... tangent function on [-pi/4,pi/4]
* __ieee754_rem_pio2 ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
* in [-pi/4 , +pi/4], and let n = k mod 4.
* We have
*
* n sin(x) cos(x) tan(x)
* ----------------------------------------------------------
* 0 S C T
* 1 C -S -1/T
* 2 -S -C T
* 3 -C S -1/T
* ----------------------------------------------------------
*
* Special cases:
* Let trig be any of sin, cos, or tan.
* trig(+-INF) is NaN, with signals;
* trig(NaN) is that NaN;
*
* Accuracy:
* TRIG(x) returns trig(x) nearly rounded
*/
double math_tan(double x)
{
double y[2],z=0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if(ix <= 0x3fe921fb) return __kernel_tan(x,z,1);
/* tan(Inf or NaN) is NaN */
else if (ix>=0x7ff00000) return x-x; /* NaN */
/* argument reduction needed */
else {
n = __ieee754_rem_pio2(x,y);
return __kernel_tan(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even
-1 -- n odd */
}
}
/* __ieee754_asin(x)
* Method :
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
* we approximate asin(x) on [0,0.5] by
* asin(x) = x + x*x^2*R(x^2)
* where
* R(x^2) is a rational approximation of (asin(x)-x)/x^3
* and its remez error is bounded by
* |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
*
* For x in [0.5,1]
* asin(x) = pi/2-2*asin(sqrt((1-x)/2))
* Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
* then for x>0.98
* asin(x) = pi/2 - 2*(s+s*z*R(z))
* = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
* For x<=0.98, let pio4_hi = pio2_hi/2, then
* f = hi part of s;
* c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
* and
* asin(x) = pi/2 - 2*(s+s*z*R(z))
* = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
* = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
*
* Special cases:
* if x is NaN, return x itself;
* if |x|>1, return NaN with invalid signal.
*
*/
double __ieee754_asin(double x)
{
double t,w,p,q,c,r,s;
int hx,ix;
hx = __HI(x);
ix = hx&0x7fffffff;
if(ix>= 0x3ff00000) { /* |x|>= 1 */
if(((ix-0x3ff00000)|__LO(x))==0)
/* asin(1)=+-pi/2 with inexact */
return x*pio2_hi+x*pio2_lo;
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix<0x3fe00000) { /* |x|<0.5 */
if(ix<0x3e400000) { /* if |x| < 2**-27 */
if(huge+x>one) return x;/* return x with inexact if x!=0*/
} else
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x+x*w;
}
/* 1> |x|>= 0.5 */
w = one-fabs(x);
t = w*0.5;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
s = sqrt(t);
if(ix>=0x3FEF3333) { /* if |x| > 0.975 */
w = p/q;
t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
} else {
w = s;
__LO(w) = 0;
c = (t-w*w)/(s+w);
r = p/q;
p = 2.0*s*r-(pio2_lo-2.0*c);
q = pio4_hi-2.0*w;
t = pio4_hi-(p-q);
}
if(hx>0) return t; else return -t;
}
/*
* wrapper asin(x)
*/
double math_asin(double x) /* wrapper asin */
{
#ifdef _IEEE_LIBM
return __ieee754_asin(x);
#else
double z;
z = __ieee754_asin(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(fabs(x)>1.0) {
return __kernel_standard(x,x,2); /* asin(|x|>1) */
} else
return z;
#endif
}
/* __ieee754_acos(x)
* Method :
* acos(x) = pi/2 - asin(x)
* acos(-x) = pi/2 + asin(x)
* For |x|<=0.5
* acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
* For x>0.5
* acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
* = 2asin(sqrt((1-x)/2))
* = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
* = 2f + (2c + 2s*z*R(z))
* where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
* for f so that f+c ~ sqrt(z).
* For x<-0.5
* acos(x) = pi - 2asin(sqrt((1-|x|)/2))
* = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
*
* Special cases:
* if x is NaN, return x itself;
* if |x|>1, return NaN with invalid signal.
*
* Function needed: sqrt
*/
double __ieee754_acos(double x)
{
double z,p,q,r,w,s,c,df;
int hx,ix;
hx = __HI(x);
ix = hx&0x7fffffff;
if(ix>=0x3ff00000) { /* |x| >= 1 */
if(((ix-0x3ff00000)|__LO(x))==0) { /* |x|==1 */
if(hx>0) return 0.0; /* acos(1) = 0 */
else return pi+2.0*pio2_lo; /* acos(-1)= pi */
}
return (x-x)/(x-x); /* acos(|x|>1) is NaN */
}
if(ix<0x3fe00000) { /* |x| < 0.5 */
if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
z = x*x;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx<0) { /* x < -0.5 */
z = (one+x)*0.5;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
s = sqrt(z);
r = p/q;
w = r*s-pio2_lo;
return pi - 2.0*(s+w);
} else { /* x > 0.5 */
z = (one-x)*0.5;
s = sqrt(z);
df = s;
__LO(df) = 0;
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
w = r*s+c;
return 2.0*(df+w);
}
}
/*
* wrap_acos(x)
*/
double math_acos(double x) /* wrapper acos */
{
#ifdef _IEEE_LIBM
return __ieee754_acos(x);
#else
double z;
z = __ieee754_acos(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(fabs(x)>1.0) {
return __kernel_standard(x,x,1); /* acos(|x|>1) */
} else
return z;
#endif
}
/* atan(x)
* Method
* 1. Reduce x to positive by atan(x) = -atan(-x).
* 2. According to the integer k=4t+0.25 chopped, t=x, the argument
* is further reduced to one of the following intervals and the
* arctangent of t is evaluated by the corresponding formula:
*
* [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
* [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
* [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
* [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
* [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
static const double atanhi[] = {
4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
};
static const double atanlo[] = {
2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
};
static const double aT[] = {
3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
-1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
-1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
-7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
-5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
-3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
};
double math_atan(double x)
{
double w,s1,s2,z;
int ix,hx,id;
hx = __HI(x);
ix = hx&0x7fffffff;
if(ix>=0x44100000) { /* if |x| >= 2^66 */
if(ix>0x7ff00000||
(ix==0x7ff00000&&(__LO(x)!=0)))
return x+x; /* NaN */
if(hx>0) return atanhi[3]+atanlo[3];
else return -atanhi[3]-atanlo[3];
} if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e200000) { /* |x| < 2^-29 */
if(huge+x>one) return x; /* raise inexact */
}
id = -1;
} else {
x = fabs(x);
if (ix < 0x3ff30000) { /* |x| < 1.1875 */
if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
id = 0; x = (2.0*x-one)/(2.0+x);
} else { /* 11/16<=|x|< 19/16 */
id = 1; x = (x-one)/(x+one);
}
} else {
if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2; x = (x-1.5)/(one+1.5*x);
} else { /* 2.4375 <= |x| < 2^66 */
id = 3; x = -1.0/x;
}
}}
/* end of argument reduction */
z = x*x;
w = z*z;
/* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
if (id<0) return x - x*(s1+s2);
else {
z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
return (hx<0)? -z:z;
}
}
/* __ieee754_sinh(x)
* Method :
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
* 1. Replace x by |x| (sinh(-x) = -sinh(x)).
* 2.
* E + E/(E+1)
* 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
* 2
*
* 22 <= x <= lnovft : sinh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : sinh(x) := x*shuge (overflow)
*
* Special cases:
* sinh(x) is |x| if x is +INF, -INF, or NaN.
* only sinh(0)=0 is exact for finite x.
*/
double __ieee754_sinh(double x)
{
double t,w,h;
int ix,jx;
unsigned lx;
/* High word of |x|. */
jx = __HI(x);
ix = jx&0x7fffffff;
/* x is INF or NaN */
if(ix>=0x7ff00000) 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 (ix<0x3e300000) /* |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 [log(maxdouble), overflowthresold] */
lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
if ( ix<0x408633CE || ((ix==0x408633ce) &&(lx<=(unsigned)0x8fb9f87d))) {
w = __ieee754_exp(0.5*fabs(x));
t = h*w;
return t*w;
}
/* |x| > overflowthresold, sinh(x) overflow */
return x*shuge;
}
/*
* wrapper sinh(x)
*/
double math_sinh(double x) /* wrapper sinh */
{
#ifdef _IEEE_LIBM
return __ieee754_sinh(x);
#else
double z;
z = __ieee754_sinh(x);
if(_LIB_VERSION == _IEEE_) return z;
if(!finite(z)&&finite(x)) {
return __kernel_standard(x,x,25); /* sinh overflow */
} else
return z;
#endif
}
/* __ieee754_cosh(x)
* Method :
* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (cosh(x) = cosh(-x)).
* 2.
* [ exp(x) - 1 ]^2
* 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
* 2*exp(x)
*
* exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : cosh(x) := -------------------
* 2
* 22 <= x <= lnovft : cosh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : cosh(x) := huge*huge (overflow)
*
* Special cases:
* cosh(x) is |x| if x is +INF, -INF, or NaN.
* only cosh(0)=1 is exact for finite x.
*/
double __ieee754_cosh(double x)
{
double t,w;
int ix;
unsigned lx;
/* High word of |x|. */
ix = __HI(x);
ix &= 0x7fffffff;
/* x is INF or NaN */
if(ix>=0x7ff00000) return x*x;
/* |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; */
if (ix < 0x40360000) {
t = __ieee754_exp(fabs(x));
return half*t+half/t;
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
if (ix<0x408633CE ||
((ix==0x408633ce)&&(lx<=(unsigned)0x8fb9f87d))) {
w = __ieee754_exp(half*fabs(x));
t = half*w;
return t*w;
}
/* |x| > overflowthresold, cosh(x) overflow */
return huge*huge;
}
/*
* wrapper cosh(x)
*/
double math_cosh(double x) /* wrapper cosh */
{
#ifdef _IEEE_LIBM
return __ieee754_cosh(x);
#else
double z;
z = __ieee754_cosh(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(fabs(x)>7.10475860073943863426e+02) {
return __kernel_standard(x,x,5); /* cosh overflow */
} else
return z;
#endif
}
/* Tanh(x)
* Return the Hyperbolic Tangent of x
*
* Method :
* x -x
* e - e
* 0. tanh(x) is defined to be -----------
* x -x
* e + e
* 1. reduce x to non-negative by tanh(-x) = -tanh(x).
* 2. 0 <= x <= 2**-55 : tanh(x) := x*(one+x)
* -t
* 2**-55 < x <= 1 : tanh(x) := -----; t = expm1(-2x)
* t + 2
* 2
* 1 <= x <= 22.0 : tanh(x) := 1- ----- ; t=expm1(2x)
* t + 2
* 22.0 < x <= INF : tanh(x) := 1.
*
* Special cases:
* tanh(NaN) is NaN;
* only tanh(0)=0 is exact for finite argument.
*/
double math_tanh(double x)
{
double t,z;
int jx,ix;
/* High word of |x|. */
jx = __HI(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| < 22 */
if (ix < 0x40360000) { /* |x|<22 */
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 */
}
return (jx>=0)? z: -z;
}
/* asinh(x)
* Method :
* Based on
* asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
* we have
* asinh(x) := x if 1+x*x=1,
* := sign(x)*(log(x)+ln2)) for large |x|, else
* := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
double math_asinh(double x)
{
double t,w;
int hx,ix;
hx = __HI(x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */
if(ix< 0x3e300000) { /* |x|<2**-28 */
if(huge+x>one) return x; /* return x inexact except 0 */
}
if(ix>0x41b00000) { /* |x| > 2**28 */
w = __ieee754_log(fabs(x))+ln2;
} else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabs(x);
w = __ieee754_log(2.0*t+one/(sqrt(x*x+one)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1p(fabs(x)+t/(one+sqrt(one+t)));
}
if(hx>0) return w; else return -w;
}
/* __ieee754_acosh(x)
* Method :
* Based on
* acosh(x) = log [ x + sqrt(x*x-1) ]
* we have
* acosh(x) := log(x)+ln2, if x is large; else
* acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
* acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
*
* Special cases:
* acosh(x) is NaN with signal if x<1.
* acosh(NaN) is NaN without signal.
*/
double __ieee754_acosh(double x)
{
double t;
int hx;
hx = __HI(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)|__LO(x))==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+sqrt(t-one)));
} else { /* 1<x<2 */
t = x-one;
return log1p(t+sqrt(2.0*t+t*t));
}
}
/*
* wrapper acosh(x)
*/
double math_acosh(double x) /* wrapper acosh */
{
#ifdef _IEEE_LIBM
return __ieee754_acosh(x);
#else
double z;
z = __ieee754_acosh(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(x<1.0) {
return __kernel_standard(x,x,29); /* acosh(x<1) */
} else
return z;
#endif
}
/* __ieee754_atanh(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
*
* Special cases:
* atanh(x) is NaN if |x| > 1 with signal;
* atanh(NaN) is that NaN with no signal;
* atanh(+-1) is +-INF with signal.
*
*/
double __ieee754_atanh(double x)
{
double t;
int hx,ix;
unsigned lx;
hx = __HI(x); /* high word */
lx = __LO(x); /* low word */
ix = hx&0x7fffffff;
if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
return (x-x)/(x-x);
if(ix==0x3ff00000)
return x/zero;
if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */
__HI(x) = ix; /* x <- |x| */
if(ix<0x3fe00000) { /* x < 0.5 */
t = x+x;
t = 0.5*log1p(t+t*x/(one-x));
} else
t = 0.5*log1p((x+x)/(one-x));
if(hx>=0) return t; else return -t;
}
/*
* wrapper atanh(x)
*/
double math_atanh(double x) /* wrapper atanh */
{
#ifdef _IEEE_LIBM
return __ieee754_atanh(x);
#else
double z,y;
z = __ieee754_atanh(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
y = fabs(x);
if(y>=1.0) {
if(y>1.0)
return __kernel_standard(x,x,30); /* atanh(|x|>1) */
else
return __kernel_standard(x,x,31); /* atanh(|x|==1) */
} else
return z;
#endif
}
/*
* wrapper exp(x)
*/
double math_exp(double x) /* wrapper exp */
{
#ifdef _IEEE_LIBM
return __ieee754_exp(x);
#else
double z;
z = __ieee754_exp(x);
if(_LIB_VERSION == _IEEE_) return z;
if(finite(x)) {
if(x>o_threshold)
return __kernel_standard(x,x,6); /* exp overflow */
else if(x<u_threshold)
return __kernel_standard(x,x,7); /* exp underflow */
}
return z;
#endif
}
/*
* fabs(x) returns the absolute value of x.
*/
double fabs(double x)
{
__HI(x) &= 0x7fffffff;
return x;
}
/*
* wrapper log(x)
*/
double math_log(double x) /* wrapper log */
{
#ifdef _IEEE_LIBM
return __ieee754_log(x);
#else
double z;
z = __ieee754_log(x);
if(_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0) return z;
if(x==0.0)
return __kernel_standard(x,x,16); /* log(0) */
else
return __kernel_standard(x,x,17); /* log(x<0) */
#endif
}
/* __ieee754_log10(x)
* Return the base 10 logarithm of x
*
* Method :
* Let log10_2hi = leading 40 bits of log10(2) and
* log10_2lo = log10(2) - log10_2hi,
* ivln10 = 1/log(10) rounded.
* Then
* n = ilogb(x),
* if(n<0) n = n+1;
* x = scalbn(x,-n);
* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
*
* Note 1:
* To guarantee log10(10**n)=n, where 10**n is normal, the rounding
* mode must set to Round-to-Nearest.
* Note 2:
* [1/log(10)] rounded to 53 bits has error .198 ulps;
* log10 is monotonic at all binary break points.
*
* Special cases:
* log10(x) is NaN with signal if x < 0;
* log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
* log10(NaN) is that NaN with no signal;
* log10(10**N) = N for N=0,1,...,22.
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
double __ieee754_log10(double x)
{
double y,z;
int i,k,hx;
unsigned lx;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
hx = __HI(x); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
i = ((unsigned)k&0x80000000)>>31;
hx = (hx&0x000fffff)|((0x3ff-i)<<20);
y = (double)(k+i);
__HI(x) = hx;
z = y*log10_2lo + ivln10*__ieee754_log(x);
return z+y*log10_2hi;
}
/*
* wrapper log10(X)
*/
double math_log10(double x) /* wrapper log10 */
{
#ifdef _IEEE_LIBM
return __ieee754_log10(x);
#else
double z;
z = __ieee754_log10(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(x<=0.0) {
if(x==0.0)
return __kernel_standard(x,x,18); /* log10(0) */
else
return __kernel_standard(x,x,19); /* log10(x<0) */
} else
return z;
#endif
}
/* __ieee754_pow(x,y) return x**y
*
* n
* Method: Let x = 2 * (1+f)
* 1. Compute and return log2(x) in two pieces:
* log2(x) = w1 + w2,
* where w1 has 53-24 = 29 bit trailing zeros.
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'|<=0.5.
* 3. Return x**y = 2**n*exp(y'*log2)
*
* Special cases:
* 1. (anything) ** 0 is 1
* 2. (anything) ** 1 is itself
* 3. (anything) ** NAN is NAN
* 4. NAN ** (anything except 0) is NAN
* 5. +-(|x| > 1) ** +INF is +INF
* 6. +-(|x| > 1) ** -INF is +0
* 7. +-(|x| < 1) ** +INF is +0
* 8. +-(|x| < 1) ** -INF is +INF
* 9. +-1 ** +-INF is NAN
* 10. +0 ** (+anything except 0, NAN) is +0
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
* 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
* 15. +INF ** (+anything except 0,NAN) is +INF
* 16. +INF ** (-anything except 0,NAN) is +0
* 17. -INF ** (anything) = -0 ** (-anything)
* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
* 19. (-anything except 0 and inf) ** (non-integer) is NAN
*
* Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular
* pow(integer,integer)
* always returns the correct integer provided it is
* representable.
*
* Constants :
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
double __ieee754_pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;
int i0,i1,i,j,k,yisint,n;
int hx,hy,ix,iy;
unsigned lx,ly;
i0 = ((*(int*)&one)>>29)^1; i1=1-i0;
hx = __HI(x); lx = __LO(x);
hy = __HI(y); ly = __LO(y);
ix = hx&0x7fffffff; iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<0) {
if(iy>=0x43400000) yisint = 2; /* even integer y */
else if(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */
if(k>20) {
j = ly>>(52-k);
if((j<<(52-k))==ly) yisint = 2-(j&1);
} else if(ly==0) {
j = iy>>(20-k);
if((j<<(20-k))==iy) yisint = 2-(j&1);
}
}
}
/* special value of y */
if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if(((ix-0x3ff00000)|lx)==0)
return y - y; /* inf**+-1 is NaN */
else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */
return (hy<0)?-y: zero;
}
if(iy==0x3ff00000) { /* y is +-1 */
if(hy<0) return one/x; else return x;
}
if(hy==0x40000000) return x*x; /* y is 2 */
if(hy==0x3fe00000) { /* y is 0.5 */
if(hx>=0) /* x >= +0 */
return sqrt(x);
}
}
ax = fabs(x);
/* special value of x */
if(lx==0) {
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
z = ax; /*x is +-0,+-inf,+-1*/
if(hy<0) z = one/z; /* z = (1/|x|) */
if(hx<0) {
if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} else if(yisint==1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
n = (hx>>31)+1;
/* (x<0)**(non-int) is NaN */
if((n|yisint)==0) return (x-x)/(x-x);
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
/* |y| is huge */
if(iy>0x41e00000) { /* if |y| > 2**31 */
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
}
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;
/* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax-one; /* t has 20 trailing zeros */
w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
__LO(t1) = 0;
t2 = v-(t1-u);
} else {
double ss,s2,s_h,s_l,t_h,t_l;
n = 0;
/* take care subnormal number */
if(ix<0x00100000)
{ax *= two53; n -= 53; ix = __HI(ax); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
ix = j|0x3ff00000; /* normalize ix */
if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
else {k=0;n+=1;ix -= 0x00100000;}
__HI(ax) = ix;
/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
ss = u*v;
s_h = ss;
__LO(s_h) = 0;
/* t_h=ax+bp[k] High */
t_h = zero;
__HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18);
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
s2 = ss*ss;
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+ss);
s2 = s_h*s_h;
t_h = 3.0+s2+r;
__LO(t_h) = 0;
t_l = r-((t_h-3.0)-s2);
/* u+v = ss*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*ss;
/* 2/(3log2)*(ss+...) */
p_h = u+v;
__LO(p_h) = 0;
p_l = v-(p_h-u);
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (double)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
__LO(t1) = 0;
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
__LO(y1) = 0;
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
j = __HI(z);
i = __LO(z);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
else {
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
}
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return s*tiny*tiny; /* underflow */
else {
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
}
}
/*
* compute 2**(p_h+p_l)
*/
i = j&0x7fffffff;
k = (i>>20)-0x3ff;
n = 0;
if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
n = j+(0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
t = zero;
__HI(t) = (n&~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
__LO(t) = 0;
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
w = v-(z-u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two)-(w+z*w);
z = one-(r-z);
j = __HI(z);
j += (n<<20);
if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
else __HI(z) += (n<<20);
return s*z;
}
/*
* wrapper pow(x,y) return x**y
*/
double math_pow(double x, double y) /* wrapper pow */
{
#ifdef _IEEE_LIBM
return __ieee754_pow(x,y);
#else
double z;
z=__ieee754_pow(x,y);
if(_LIB_VERSION == _IEEE_|| isnan(y)) return z;
if(isnan(x)) {
if(y==0.0)
return __kernel_standard(x,y,42); /* pow(NaN,0.0) */
else
return z;
}
if(x==0.0){
if(y==0.0)
return __kernel_standard(x,y,20); /* pow(0.0,0.0) */
if(finite(y)&&y<0.0)
return __kernel_standard(x,y,23); /* pow(0.0,negative) */
return z;
}
if(!finite(z)) {
if(finite(x)&&finite(y)) {
if(isnan(z))
return __kernel_standard(x,y,24); /* pow neg**non-int */
else
return __kernel_standard(x,y,21); /* pow overflow */
}
}
if(z==0.0&&finite(x)&&finite(y))
return __kernel_standard(x,y,22); /* pow underflow */
return z;
#endif
}
/* __ieee754_sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
* | Use the hardware sqrt if you have one |
* ------------------------------------------
* Method:
* Bit by bit method using integer arithmetic. (Slow, but portable)
* 1. Normalization
* Scale x to y in [1,4) with even powers of 2:
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
* sqrt(x) = 2^k * sqrt(y)
* 2. Bit by bit computation
* Let q = sqrt(y) truncated to i bit after binary point (q = 1),
* i 0
* i+1 2
* s = 2*q , and y = 2 * ( y - q ). (1)
* i i i i
*
* To compute q from q , one checks whether
* i+1 i
*
* -(i+1) 2
* (q + 2 ) <= y. (2)
* i
* -(i+1)
* If (2) is false, then q = q ; otherwise q = q + 2 .
* i+1 i i+1 i
*
* With some algebric manipulation, it is not difficult to see
* that (2) is equivalent to
* -(i+1)
* s + 2 <= y (3)
* i i
*
* The advantage of (3) is that s and y can be computed by
* i i
* the following recurrence formula:
* if (3) is false
*
* s = s , y = y ; (4)
* i+1 i i+1 i
*
* otherwise,
* -i -(i+1)
* s = s + 2 , y = y - s - 2 (5)
* i+1 i i+1 i i
*
* One may easily use induction to prove (4) and (5).
* Note. Since the left hand side of (3) contain only i+2 bits,
* it does not necessary to do a full (53-bit) comparison
* in (3).
* 3. Final rounding
* After generating the 53 bits result, we compute one more bit.
* Together with the remainder, we can decide whether the
* result is exact, bigger than 1/2ulp, or less than 1/2ulp
* (it will never equal to 1/2ulp).
* The rounding mode can be detected by checking whether
* huge + tiny is equal to huge, and whether huge - tiny is
* equal to huge for some floating point number "huge" and "tiny".
*
* Special cases:
* sqrt(+-0) = +-0 ... exact
* sqrt(inf) = inf
* sqrt(-ve) = NaN ... with invalid signal
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
*
* Other methods : see the appended file at the end of the program below.
*---------------
*/
double __ieee754_sqrt(double x)
{
double z;
int sign = (int)0x80000000;
unsigned r,t1,s1,ix1,q1;
int ix0,s0,q,m,t,i;
ix0 = __HI(x); /* high word of x */
ix1 = __LO(x); /* low word of x */
/* take care of Inf and NaN */
if((ix0&0x7ff00000)==0x7ff00000) {
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
sqrt(-inf)=sNaN */
}
/* take care of zero */
if(ix0<=0) {
if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
else if(ix0<0)
return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
}
/* normalize x */
m = (ix0>>20);
if(m==0) { /* subnormal x */
while(ix0==0) {
m -= 21;
ix0 |= (ix1>>11); ix1 <<= 21;
}
for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
m -= i-1;
ix0 |= (ix1>>(32-i));
ix1 <<= i;
}
m -= 1023; /* unbias exponent */
ix0 = (ix0&0x000fffff)|0x00100000;
if(m&1){ /* odd m, double x to make it even */
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
r = 0x00200000; /* r = moving bit from right to left */
while(r!=0) {
t = s0+r;
if(t<=ix0) {
s0 = t+r;
ix0 -= t;
q += r;
}
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
r>>=1;
}
r = sign;
while(r!=0) {
t1 = s1+r;
t = s0;
if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
s1 = t1+r;
if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
ix0 -= t;
if (ix1 < t1) ix0 -= 1;
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1&sign)>>31);
ix1 += ix1;
r>>=1;
}
/* use floating add to find out rounding direction */
if((ix0|ix1)!=0) {
z = one-tiny; /* trigger inexact flag */
if (z>=one) {
z = one+tiny;
if (q1==(unsigned)0xffffffff) { q1=0; q += 1;}
else if (z>one) {
if (q1==(unsigned)0xfffffffe) q+=1;
q1+=2;
} else
q1 += (q1&1);
}
}
ix0 = (q>>1)+0x3fe00000;
ix1 = q1>>1;
if ((q&1)==1) ix1 |= sign;
ix0 += (m <<20);
__HI(z) = ix0;
__LO(z) = ix1;
return z;
}
/*
Other methods (use floating-point arithmetic)
-------------
(This is a copy of a drafted paper by Prof W. Kahan
and K.C. Ng, written in May, 1986)
Two algorithms are given here to implement sqrt(x)
(IEEE double precision arithmetic) in software.
Both supply sqrt(x) correctly rounded. The first algorithm (in
Section A) uses newton iterations and involves four divisions.
The second one uses reciproot iterations to avoid division, but
requires more multiplications. Both algorithms need the ability
to chop results of arithmetic operations instead of round them,
and the INEXACT flag to indicate when an arithmetic operation
is executed exactly with no roundoff error, all part of the
standard (IEEE 754-1985). The ability to perform shift, add,
subtract and logical AND operations upon 32-bit words is needed
too, though not part of the standard.
A. sqrt(x) by Newton Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
1 11 52 ...widths
------------------------------------------------------
x: |s| e | f |
------------------------------------------------------
msb lsb msb lsb ...order
------------------------ ------------------------
x0: |s| e | f1 | x1: | f2 |
------------------------ ------------------------
By performing shifts and subtracts on x0 and x1 (both regarded
as integers), we obtain an 8-bit approximation of sqrt(x) as
follows.
k := (x0>>1) + 0x1ff80000;
y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
Here k is a 32-bit integer and T1[] is an integer array containing
correction terms. Now magically the floating value of y (y's
leading 32-bit word is y0, the value of its trailing word is 0)
approximates sqrt(x) to almost 8-bit.
Value of T1:
static int T1[32]= {
0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
(2) Iterative refinement
Apply Heron's rule three times to y, we have y approximates
sqrt(x) to within 1 ulp (Unit in the Last Place):
y := (y+x/y)/2 ... almost 17 sig. bits
y := (y+x/y)/2 ... almost 35 sig. bits
y := y-(y-x/y)/2 ... within 1 ulp
Remark 1.
Another way to improve y to within 1 ulp is:
y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
2
(x-y )*y
y := y + 2* ---------- ...within 1 ulp
2
3y + x
This formula has one division fewer than the one above; however,
it requires more multiplications and additions. Also x must be
scaled in advance to avoid spurious overflow in evaluating the
expression 3y*y+x. Hence it is not recommended uless division
is slow. If division is very slow, then one should use the
reciproot algorithm given in section B.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
I := FALSE; ... reset INEXACT flag I
R := RZ; ... set rounding mode to round-toward-zero
z := x/y; ... chopped quotient, possibly inexact
If(not I) then { ... if the quotient is exact
if(z=y) {
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
} else {
z := z - ulp; ... special rounding
}
}
i := TRUE; ... sqrt(x) is inexact
If (r=RN) then z=z+ulp ... rounded-to-nearest
If (r=RP) then { ... round-toward-+inf
y = y+ulp; z=z+ulp;
}
y := y+z; ... chopped sum
y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
(4) Special cases
Square root of +inf, +-0, or NaN is itself;
Square root of a negative number is NaN with invalid signal.
B. sqrt(x) by Reciproot Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
(see section A). By performing shifs and subtracts on x0 and y0,
we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
k := 0x5fe80000 - (x0>>1);
y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
Here k is a 32-bit integer and T2[] is an integer array
containing correction terms. Now magically the floating
value of y (y's leading 32-bit word is y0, the value of
its trailing word y1 is set to zero) approximates 1/sqrt(x)
to almost 7.8-bit.
Value of T2:
static int T2[64]= {
0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
(2) Iterative refinement
Apply Reciproot iteration three times to y and multiply the
result by x to get an approximation z that matches sqrt(x)
to about 1 ulp. To be exact, we will have
-1ulp < sqrt(x)-z<1.0625ulp.
... set rounding mode to Round-to-nearest
y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
... special arrangement for better accuracy
z := x*y ... 29 bits to sqrt(x), with z*y<1
z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
(a) the term z*y in the final iteration is always less than 1;
(b) the error in the final result is biased upward so that
-1 ulp < sqrt(x) - z < 1.0625 ulp
instead of |sqrt(x)-z|<1.03125ulp.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
R := RZ; ... set rounding mode to round-toward-zero
switch(r) {
case RN: ... round-to-nearest
if(x<= z*(z-ulp)...chopped) z = z - ulp; else
if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
break;
case RZ:case RM: ... round-to-zero or round-to--inf
R:=RP; ... reset rounding mod to round-to-+inf
if(x<z*z ... rounded up) z = z - ulp; else
if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
break;
case RP: ... round-to-+inf
if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
if(x>z*z ...chopped) z = z+ulp;
break;
}
Remark 3. The above comparisons can be done in fixed point. For
example, to compare x and w=z*z chopped, it suffices to compare
x1 and w1 (the trailing parts of x and w), regarding them as
two's complement integers.
...Is z an exact square root?
To determine whether z is an exact square root of x, let z1 be the
trailing part of z, and also let x0 and x1 be the leading and
trailing parts of x.
If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
I := 1; ... Raise Inexact flag: z is not exact
else {
j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
k := z1 >> 26; ... get z's 25-th and 26-th
fraction bits
I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
}
R:= r ... restore rounded mode
return sqrt(x):=z.
If multiplication is cheaper then the foregoing red tape, the
Inexact flag can be evaluated by
I := i;
I := (z*z!=x) or I.
Note that z*z can overwrite I; this value must be sensed if it is
True.
Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
zero.
--------------------
z1: | f2 |
--------------------
bit 31 bit 0
Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
or even of logb(x) have the following relations:
-------------------------------------------------
bit 27,26 of z1 bit 1,0 of x1 logb(x)
-------------------------------------------------
00 00 odd and even
01 01 even
10 10 odd
10 00 even
11 01 even
-------------------------------------------------
(4) Special cases (see (4) of Section A).
*/
/*
* wrapper sqrt(x)
*/
double math_sqrt(double x) /* wrapper sqrt */
{
#ifdef _IEEE_LIBM
return __ieee754_sqrt(x);
#else
double z;
z = __ieee754_sqrt(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(x<0.0) {
return __kernel_standard(x,x,26); /* sqrt(negative) */
} else
return z;
#endif
}
/*
* ceil(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to ceil(x).
*/
double math_ceil(double x)
{
int i0,i1,j0;
unsigned i,j;
i0 = __HI(x);
i1 = __LO(x);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* 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 */
if(huge+x>0.0) { /* 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 = ((unsigned)(0xffffffff))>>(j0-20);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* 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);
}
}
__HI(x) = i0;
__LO(x) = i1;
return x;
}
/*
* floor(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to floor(x).
*/
double math_floor(double x)
{
int i0,i1,j0;
unsigned i,j;
i0 = __HI(x);
i1 = __LO(x);
j0 = ((i0>>20)&0x7ff)-0x3ff;
if(j0<20) {
if(j0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* 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 */
if(huge+x>0.0) { /* 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 = ((unsigned)(0xffffffff))>>(j0-20);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) { /* 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);
}
}
__HI(x) = i0;
__LO(x) = i1;
return x;
}
/*
* fabs(x) returns the absolute value of x.
*/
double math_fabs(double x)
{
__HI(x) &= 0x7fffffff;
return x;
}
#endif /* NEED_MATH_LIBRARY */
Reference in a new issue View git blame Copy permalink