zero, one, two, half are replaced by const literals
The policy was to use the f suffix for float consts (1.0f),
but don't use suffix for long double consts (these consts
can be exactly represented as double).
73 files changed:
#include "libm.h"
static const double
#include "libm.h"
static const double
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
w = z*z;
r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
hz = 0.5*z;
w = z*z;
r = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
hz = 0.5*z;
- w = one-hz;
- return w + (((one-w)-hz) + (z*r-x*y));
+ w = 1.0-hz;
+ return w + (((1.0-w)-hz) + (z*r-x*y));
/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
static const double
/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
static const double
C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
z = x*x;
w = z*z;
r = C2+z*C3;
z = x*x;
w = z*z;
r = C2+z*C3;
- return ((one+z*C0) + w*C1) + (w*z)*r;
+ return ((1.0+z*C0) + w*C1) + (w*z)*r;
* almost for free from the complications needed to search for the best
* higher coefficients.
*/
* almost for free from the complications needed to search for the best
* higher coefficients.
*/
-static const double one = 1.0;
-
static const long double
C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
static const long double
C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))));
hz = 0.5*z;
z = x*x;
r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))));
hz = 0.5*z;
- w = one-hz;
- return w + (((one-w)-hz) + (z*r-x*y));
+ w = 1.0-hz;
+ return w + (((1.0-w)-hz) + (z*r-x*y));
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
-zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
- while (tx[nx-1] == zero) nx--; /* skip zero term */
+ while (tx[nx-1] == 0.0) nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,1);
if (hx < 0) {
y[0] = -ty[0];
n = __rem_pio2_large(tx,ty,e0,nx,1);
if (hx < 0) {
y[0] = -ty[0];
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
/* 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++)
/* 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];
+ f[i] = j<0 ? 0.0 : (double)ipio2[j];
/* compute q[0],q[1],...q[jk] */
for (i=0; i<=jk; i++) {
/* compute q[0],q[1],...q[jk] */
for (i=0; i<=jk; i++) {
}
}
/* check if recomputation is needed */
}
}
/* check if recomputation is needed */
j = 0;
for (i=jz-1; i>=jk; i--) j |= iq[i];
if (j == 0) { /* need recomputation */
j = 0;
for (i=jz-1; i>=jk; i--) j |= iq[i];
if (j == 0) { /* need recomputation */
}
/* convert integer "bit" chunk to floating-point value */
}
/* convert integer "bit" chunk to floating-point value */
for (i=jz; i>=0; i--) {
q[i] = fw*(double)iq[i];
fw *= twon24;
for (i=jz; i>=0; i--) {
q[i] = fw*(double)iq[i];
fw *= twon24;
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
static const double
-zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
- while (tx[nx-1] == zero)
+ while (tx[nx-1] == 0.0)
nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,2);
r = (long double)ty[0] + ty[1];
nx--; /* skip zero term */
n = __rem_pio2_large(tx,ty,e0,nx,2);
r = (long double)ty[0] + ty[1];
#include "libm.h"
static const double
#include "libm.h"
static const double
-half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
if (iy == 0)
return x + v*(S1 + z*r);
else
if (iy == 0)
return x + v*(S1 + z*r);
else
- return x - ((z*(half*y - v*r) - y) - v*S1);
+ return x - ((z*(0.5*y - v*r) - y) - v*S1);
* See __cosl.c for more details about the polynomial.
*/
* See __cosl.c for more details about the polynomial.
*/
-static const double half = 0.5;
-
static const long double
S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
static const long double
S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))));
if (iy == 0)
return x+v*(S1+z*r);
r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))));
if (iy == 0)
return x+v*(S1+z*r);
- return x-((z*(half*y-v*r)-y)-v*S1);
+ return x-((z*(0.5*y-v*r)-y)-v*S1);
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */
-1.85586374855275456654e-05, /* BEF375CB, DB605373 */
2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */
-/* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */
-/* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
-/* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */
-};
-#define one T[13]
-#define pio4 T[14]
-#define pio4lo T[15]
+},
+pio4 = 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */
+pio4lo = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
double __tan(double x, double y, int iy)
{
double __tan(double x, double y, int iy)
{
#include "libm.h"
static const double
#include "libm.h"
static const double
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
static const volatile double
pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
static const volatile double
pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
return pio2_hi + pio2_lo;
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)));
+ q = 1.0+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 */
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
- q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ q = 1.0+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 */
s = sqrt(z);
r = p/q;
w = r*s-pio2_lo;
return pi - 2.0*(s+w);
} else { /* x > 0.5 */
s = sqrt(z);
df = s;
SET_LOW_WORD(df,0);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
s = sqrt(z);
df = s;
SET_LOW_WORD(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)));
+ q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
w = r*s+c;
return 2.0*(df+w);
r = p/q;
w = r*s+c;
return 2.0*(df+w);
#include "libm.h"
static const float
#include "libm.h"
static const float
-one = 1.0000000000e+00, /* 0x3F800000 */
-pi = 3.1415925026e+00, /* 0x40490fda */
+pi = 3.1415925026e+00, /* 0x40490fda */
pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
static const volatile float
pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
static const volatile float
pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*pS2));
return pio2_hi + pio2_lo;
z = x*x;
p = z*(pS0+z*(pS1+z*pS2));
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
r = p/q;
return pio2_hi - (x - (pio2_lo-x*r));
} else if (hx < 0) { /* x < -0.5 */
p = z*(pS0+z*(pS1+z*pS2));
p = z*(pS0+z*(pS1+z*pS2));
s = sqrtf(z);
r = p/q;
w = r*s-pio2_lo;
s = sqrtf(z);
r = p/q;
w = r*s-pio2_lo;
} else { /* x > 0.5 */
int32_t idf;
} else { /* x > 0.5 */
int32_t idf;
s = sqrtf(z);
df = s;
GET_FLOAT_WORD(idf,df);
SET_FLOAT_WORD(df,idf&0xfffff000);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*pS2));
s = sqrtf(z);
df = s;
GET_FLOAT_WORD(idf,df);
SET_FLOAT_WORD(df,idf&0xfffff000);
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*pS2));
r = p/q;
w = r*s+c;
return 2.0f*(df+w);
r = p/q;
w = r*s+c;
return 2.0f*(df+w);
#include "libm.h"
static const double
#include "libm.h"
static const double
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double acosh(double x)
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
double acosh(double x)
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
- return log(2.0*x - one/(x+sqrt(t-one)));
+ return log(2.0*x - 1.0/(x+sqrt(t-1.0)));
return log1p(t + sqrt(2.0*t+t*t));
}
}
return log1p(t + sqrt(2.0*t+t*t));
}
}
#include "libm.h"
static const float
#include "libm.h"
static const float
ln2 = 6.9314718246e-01; /* 0x3f317218 */
float acoshf(float x)
ln2 = 6.9314718246e-01; /* 0x3f317218 */
float acoshf(float x)
return x + x;
return logf(x) + ln2; /* acosh(huge)=log(2x) */
} else if (hx == 0x3f800000) {
return x + x;
return logf(x) + ln2; /* acosh(huge)=log(2x) */
} else if (hx == 0x3f800000) {
- return 0.0; /* acosh(1) = 0 */
+ return 0.0f; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t = x*x;
- return logf(2.0f*x - one/(x+sqrtf(t-one)));
+ return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f)));
return log1pf(t + sqrtf(2.0f*t+t*t));
}
}
return log1pf(t + sqrtf(2.0f*t+t*t));
}
}
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
long double acoshl(long double x)
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
long double acoshl(long double x)
return 0.0; /* acosh(1) = 0 */
} else if (se > 0x4000) { /* x > 2 */
t = x*x;
return 0.0; /* acosh(1) = 0 */
} else if (se > 0x4000) { /* x > 2 */
t = x*x;
- return logl(2.0*x - one/(x + sqrtl(t - one)));
+ return logl(2.0*x - 1.0/(x + sqrtl(t - 1.0)));
return log1pl(t + sqrtl(2.0*t + t*t));
}
#endif
return log1pl(t + sqrtl(2.0*t + t*t));
}
#endif
#include "__invtrigl.h"
static const long double
#include "__invtrigl.h"
static const long double
-one = 1.00000000000000000000e+00,
pi = 3.14159265358979323846264338327950280e+00L;
long double acosl(long double x)
pi = 3.14159265358979323846264338327950280e+00L;
long double acosl(long double x)
r = p / q;
return pio2_hi - (x - (pio2_lo - x * r));
} else if (expsign < 0) { /* x < -0.5 */
r = p / q;
return pio2_hi - (x - (pio2_lo - x * r));
} else if (expsign < 0) { /* x < -0.5 */
p = P(z);
q = Q(z);
s = sqrtl(z);
p = P(z);
q = Q(z);
s = sqrtl(z);
w = r * s - pio2_lo;
return pi - 2.0 * (s + w);
} else { /* x > 0.5 */
w = r * s - pio2_lo;
return pi - 2.0 * (s + w);
} else { /* x > 0.5 */
s = sqrtl(z);
u.e = s;
u.bits.manl = 0;
s = sqrtl(z);
u.e = s;
u.bits.manl = 0;
#include "libm.h"
static const double
#include "libm.h"
static const double
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
huge = 1.000e+300,
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
huge = 1.000e+300,
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3fe00000) { /* |x|<0.5 */
if (ix < 0x3e500000) { /* if |x| < 2**-26 */
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3fe00000) { /* |x|<0.5 */
if (ix < 0x3e500000) { /* if |x| < 2**-26 */
return x; /* return x with inexact if x!=0*/
}
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
return x; /* return x with inexact if x!=0*/
}
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)));
+ q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x + x*w;
}
/* 1 > |x| >= 0.5 */
w = p/q;
return x + x*w;
}
/* 1 > |x| >= 0.5 */
t = w*0.5;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
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)));
+ q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
s = sqrt(t);
if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
w = p/q;
s = sqrt(t);
if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
w = p/q;
#include "libm.h"
static const float
#include "libm.h"
static const float
-one = 1.0000000000e+00, /* 0x3F800000 */
huge = 1.000e+30,
/* coefficients for R(x^2) */
pS0 = 1.6666586697e-01,
huge = 1.000e+30,
/* coefficients for R(x^2) */
pS0 = 1.6666586697e-01,
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3f000000) { /* |x|<0.5 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
} else if (ix < 0x3f000000) { /* |x|<0.5 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
return x; /* return x with inexact if x!=0 */
}
t = x*x;
p = t*(pS0+t*(pS1+t*pS2));
return x; /* return x with inexact if x!=0 */
}
t = x*x;
p = t*(pS0+t*(pS1+t*pS2));
w = p/q;
return x + x*w;
}
/* 1 > |x| >= 0.5 */
w = p/q;
return x + x*w;
}
/* 1 > |x| >= 0.5 */
t = w*0.5f;
p = t*(pS0+t*(pS1+t*pS2));
t = w*0.5f;
p = t*(pS0+t*(pS1+t*pS2));
s = sqrt(t);
w = p/q;
t = pio2-2.0*(s+s*w);
s = sqrt(t);
w = p/q;
t = pio2-2.0*(s+s*w);
#include "libm.h"
static const double
#include "libm.h"
static const double
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300;
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300;
return x+x;
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
return x+x;
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
return x;
}
if (ix > 0x41b00000) { /* |x| > 2**28 */
w = log(fabs(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabs(x);
return x;
}
if (ix > 0x41b00000) { /* |x| > 2**28 */
w = log(fabs(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabs(x);
- w = log(2.0*t + one/(sqrt(x*x+one)+t));
+ w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
- w =log1p(fabs(x) + t/(one+sqrt(one+t)));
+ w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
#include "libm.h"
static const float
#include "libm.h"
static const float
-one = 1.0000000000e+00, /* 0x3F800000 */
ln2 = 6.9314718246e-01, /* 0x3f317218 */
huge= 1.0000000000e+30;
ln2 = 6.9314718246e-01, /* 0x3f317218 */
huge= 1.0000000000e+30;
return x+x;
if (ix < 0x31800000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
return x+x;
if (ix < 0x31800000) { /* |x| < 2**-28 */
/* return x inexact except 0 */
return x;
}
if (ix > 0x4d800000) { /* |x| > 2**28 */
w = logf(fabsf(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabsf(x);
return x;
}
if (ix > 0x4d800000) { /* |x| > 2**28 */
w = logf(fabsf(x)) + ln2;
} else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
t = fabsf(x);
- w = logf(2.0f*t + one/(sqrtf(x*x+one)+t));
+ w = logf(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
- w =log1pf(fabsf(x) + t/(one+sqrtf(one+t)));
+ w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
-one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge = 1.000000000000000000e+4900L;
ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
huge = 1.000000000000000000e+4900L;
return x + x; /* x is inf or NaN */
if (ix < 0x3fde) { /* |x| < 2**-34 */
/* return x, raise inexact if x != 0 */
return x + x; /* x is inf or NaN */
if (ix < 0x3fde) { /* |x| < 2**-34 */
/* return x, raise inexact if x != 0 */
return x;
}
if (ix > 0x4020) { /* |x| > 2**34 */
w = logl(fabsl(x)) + ln2;
} else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x);
return x;
}
if (ix > 0x4020) { /* |x| > 2**34 */
w = logl(fabsl(x)) + ln2;
} else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */
t = fabsl(x);
- w = logl(2.0*t + one/(sqrtl(x*x + one) + t));
+ w = logl(2.0*t + 1.0/(sqrtl(x*x + 1.0) + t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
- w =log1pl(fabsl(x) + t/(one + sqrtl(one + t)));
+ w =log1pl(fabsl(x) + t/(1.0 + sqrtl(1.0 + t)));
}
if (hx & 0x8000)
return -w;
}
if (hx & 0x8000)
return -w;
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
-static const long double
-one = 1.00000000000000000000e+00,
-huge = 1.000e+300;
+static const long double huge = 1.000e+300;
long double asinl(long double x)
{
long double asinl(long double x)
{
} else if (expt < BIAS-1) { /* |x|<0.5 */
if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
/* return x with inexact if x!=0 */
} else if (expt < BIAS-1) { /* |x|<0.5 */
if (expt < ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
/* return x with inexact if x!=0 */
return x + x*w;
}
/* 1 > |x| >= 0.5 */
return x + x*w;
}
/* 1 > |x| >= 0.5 */
t = w*0.5;
p = P(t);
q = Q(t);
t = w*0.5;
p = P(t);
q = Q(t);
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
};
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
};
-static const double
-one = 1.0,
-huge = 1.0e300;
+static const double huge = 1.0e300;
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact */
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e400000) { /* |x| < 2^-27 */
/* raise inexact */
if (ix < 0x3ff30000) { /* |x| < 1.1875 */
if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
id = 0;
if (ix < 0x3ff30000) { /* |x| < 1.1875 */
if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
id = 0;
- x = (2.0*x-one)/(2.0+x);
+ x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
}
} else {
if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2;
}
} else {
if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2;
- x = (x-1.5)/(one+1.5*x);
+ x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^66 */
id = 3;
x = -1.0/x;
} else { /* 2.4375 <= |x| < 2^66 */
id = 3;
x = -1.0/x;
static const volatile double
tiny = 1.0e-300;
static const double
static const volatile double
tiny = 1.0e-300;
static const double
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
- case 0: return zero; /* atan(+...,+INF) */
- case 1: return -zero; /* atan(-...,+INF) */
+ case 0: return 0.0; /* atan(+...,+INF) */
+ case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}
static const volatile float
tiny = 1.0e-30;
static const float
static const volatile float
tiny = 1.0e-30;
static const float
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
pi = 3.1415927410e+00; /* 0x40490fdb */
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
pi = 3.1415927410e+00; /* 0x40490fdb */
- case 0: return zero; /* atan(+...,+INF) */
- case 1: return -zero; /* atan(-...,+INF) */
+ case 0: return 0.0f; /* atan(+...,+INF) */
+ case 1: return -0.0f; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}
static const volatile long double
tiny = 1.0e-300;
static const long double
static const volatile long double
tiny = 1.0e-300;
static const long double
pi = 3.14159265358979323846264338327950280e+00L;
long double atan2l(long double y, long double x)
pi = 3.14159265358979323846264338327950280e+00L;
long double atan2l(long double y, long double x)
- case 0: return zero; /* atan(+...,+INF) */
- case 1: return -zero; /* atan(-...,+INF) */
+ case 0: return 0.0; /* atan(+...,+INF) */
+ case 1: return -0.0; /* atan(-...,+INF) */
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}
case 2: return pi+tiny; /* atan(+...,-INF) */
case 3: return -pi-tiny; /* atan(-...,-INF) */
}
-static const float
-one = 1.0,
-huge = 1.0e30;
+static const float huge = 1.0e30;
if (ix < 0x3ee00000) { /* |x| < 0.4375 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* raise inexact */
if (ix < 0x3ee00000) { /* |x| < 0.4375 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* raise inexact */
if (ix < 0x3f980000) { /* |x| < 1.1875 */
if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
id = 0;
if (ix < 0x3f980000) { /* |x| < 1.1875 */
if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
id = 0;
- x = (2.0f*x - one)/(2.0f + x);
+ x = (2.0f*x - 1.0f)/(2.0f + x);
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
- x = (x - one)/(x + one);
+ x = (x - 1.0f)/(x + 1.0f);
}
} else {
if (ix < 0x401c0000) { /* |x| < 2.4375 */
id = 2;
}
} else {
if (ix < 0x401c0000) { /* |x| < 2.4375 */
id = 2;
- x = (x - 1.5f)/(one + 1.5f*x);
+ x = (x - 1.5f)/(1.0f + 1.5f*x);
} else { /* 2.4375 <= |x| < 2**26 */
id = 3;
x = -1.0f/x;
} else { /* 2.4375 <= |x| < 2**26 */
id = 3;
x = -1.0f/x;
-static const double one = 1.0, huge = 1e300;
-static const double zero = 0.0;
+static const double huge = 1e300;
if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3ff00000)
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) /* x < 2**-28 */
+ return x/0.0;
+ if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */
return x;
SET_HIGH_WORD(x, ix);
if (ix < 0x3fe00000) { /* x < 0.5 */
t = x+x;
return x;
SET_HIGH_WORD(x, ix);
if (ix < 0x3fe00000) { /* x < 0.5 */
t = x+x;
- t = 0.5*log1p(t + t*x/(one-x));
+ t = 0.5*log1p(t + t*x/(1.0-x));
- t = 0.5*log1p((x+x)/(one-x));
+ t = 0.5*log1p((x+x)/(1.0-x));
if (hx >= 0)
return t;
return -t;
if (hx >= 0)
return t;
return -t;
-static const float one = 1.0, huge = 1e30;
-static const float zero = 0.0;
+static const float huge = 1e30;
if (ix > 0x3f800000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3f800000)
if (ix > 0x3f800000) /* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3f800000)
- return x/zero;
- if (ix < 0x31800000 && huge+x > zero) /* x < 2**-28 */
+ return x/0.0f;
+ if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */
return x;
SET_FLOAT_WORD(x, ix);
if (ix < 0x3f000000) { /* x < 0.5 */
t = x+x;
return x;
SET_FLOAT_WORD(x, ix);
if (ix < 0x3f000000) { /* x < 0.5 */
t = x+x;
- t = 0.5f*log1pf(t + t*x/(one-x));
+ t = 0.5f*log1pf(t + t*x/(1.0f-x));
- t = 0.5f*log1pf((x+x)/(one-x));
+ t = 0.5f*log1pf((x+x)/(1.0f-x));
if (hx >= 0)
return t;
return -t;
if (hx >= 0)
return t;
return -t;
return atanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
return atanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double zero = 0.0, one = 1.0, huge = 1e4900L;
+static const long double huge = 1e4900L;
long double atanhl(long double x)
{
long double atanhl(long double x)
{
/* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3fff)
/* |x| > 1 */
return (x-x)/(x-x);
if (ix == 0x3fff)
- return x/zero;
- if (ix < 0x3fe3 && huge+x > zero) /* x < 2**-28 */
+ return x/0.0;
+ if (ix < 0x3fe3 && huge+x > 0.0) /* x < 2**-28 */
return x;
SET_LDOUBLE_EXP(x, ix);
if (ix < 0x3ffe) { /* x < 0.5 */
t = x + x;
return x;
SET_LDOUBLE_EXP(x, ix);
if (ix < 0x3ffe) { /* x < 0.5 */
t = x + x;
- t = 0.5*log1pl(t + t*x/(one-x));
+ t = 0.5*log1pl(t + t*x/(1.0 - x));
- t = 0.5*log1pl((x + x)/(one - x));
+ t = 0.5*log1pl((x + x)/(1.0 - x));
if (se <= 0x7fff)
return t;
return -t;
if (se <= 0x7fff)
return t;
return -t;
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
-static const long double
-one = 1.0,
-huge = 1.0e300;
+static const long double huge = 1.0e300;
long double atanl(long double x)
{
long double atanl(long double x)
{
if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
/* raise inexact */
if (expman < ((BIAS - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
if (expt < ATAN_LINEAR) { /* if |x| is small, atanl(x)~=x */
/* raise inexact */
if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
id = 0;
if (expman < (BIAS << 8) + 0x30) { /* |x| < 1.1875 */
if (expman < ((BIAS - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
id = 0;
- x = (2.0*x-one)/(2.0+x);
+ x = (2.0*x-1.0)/(2.0+x);
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
} else { /* 11/16 <= |x| < 19/16 */
id = 1;
}
} else {
if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
id = 2;
}
} else {
if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
id = 2;
- x = (x-1.5)/(one+1.5*x);
+ x = (x-1.5)/(1.0+1.5*x);
} else { /* 2.4375 <= |x| < 2^ATAN_CONST */
id = 3;
x = -1.0/x;
} else { /* 2.4375 <= |x| < 2^ATAN_CONST */
id = 3;
x = -1.0/x;
-static const double one = 1.0, half = 0.5, huge = 1.0e300;
+static const double huge = 1.0e300;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43) {
t = expm1(fabs(x));
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43) {
t = expm1(fabs(x));
if (ix < 0x3c800000)
return w; /* cosh(tiny) = 1 */
if (ix < 0x3c800000)
return w; /* cosh(tiny) = 1 */
- return one + (t*t)/(w+w);
+ return 1.0 + (t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x40360000) {
t = exp(fabs(x));
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x40360000) {
t = exp(fabs(x));
- return half*t + half/t;
- /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+ /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
- return half*exp(fabs(x));
+ return 0.5*exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633CE)
/* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633CE)
-static const float one = 1.0, half = 0.5, huge = 1.0e30;
+static const float huge = 1.0e30;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3eb17218) {
t = expm1f(fabsf(x));
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3eb17218) {
t = expm1f(fabsf(x));
- return one; /* cosh(tiny) = 1 */
- return one + (t*t)/(w+w);
+ return 1.0f; /* cosh(tiny) = 1 */
+ return 1.0f + (t*t)/(w+w);
}
/* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x41100000) {
t = expf(fabsf(x));
}
/* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x41100000) {
t = expf(fabsf(x));
- return half*t + half/t;
+ return 0.5f*t + 0.5f/t;
- /* |x| in [9, log(maxfloat)] return half*exp(|x|) */
+ /* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
- return half*expf(fabsf(x));
+ return 0.5f*expf(fabsf(x));
/* |x| in [log(maxfloat), overflowthresold] */
if (ix <= 0x42b2d4fc)
/* |x| in [log(maxfloat), overflowthresold] */
if (ix <= 0x42b2d4fc)
return cosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
return cosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double one = 1.0, half = 0.5, huge = 1.0e4900L;
+static const long double huge = 1.0e4900L;
long double coshl(long double x)
{
long double coshl(long double x)
{
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
t = expm1l(fabsl(x));
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
if (ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
t = expm1l(fabsl(x));
if (ex < 0x3fbc) return w; /* cosh(tiny) = 1 */
if (ex < 0x3fbc) return w; /* cosh(tiny) = 1 */
- return one+(t*t)/(w+w);
+ return 1.0+(t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
t = expl(fabsl(x));
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
t = expl(fabsl(x));
- return half*t + half/t;
- /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
+ /* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
- return half*expl(fabsl(x));
+ return 0.5*expl(fabsl(x));
/* |x| in [log(maxdouble), log(2*maxdouble)) */
if (ex == 0x400c && (mx < 0xb174ddc0u ||
(mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
{
/* |x| in [log(maxdouble), log(2*maxdouble)) */
if (ex == 0x400c && (mx < 0xb174ddc0u ||
(mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
{
- w = expl(half*fabsl(x));
- t = half*w;
+ w = expl(0.5*fabsl(x));
+ t = 0.5*w;
static const double
tiny = 1e-300,
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 */
/*
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
if (ix >= 0x7ff00000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
if (ix >= 0x7ff00000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
- return (double)(1-i) + one/x;
+ return (double)(1-i) + 1.0/x;
}
if (ix < 0x3feb0000) { /* |x|<0.84375 */
if (ix < 0x3e300000) { /* |x|<2**-28 */
}
if (ix < 0x3feb0000) { /* |x|<0.84375 */
if (ix < 0x3e300000) { /* |x|<2**-28 */
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+ s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
y = r/s;
return x + x*y;
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+ Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0)
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf > |x| >= 6 */
if (hx >= 0)
- return one-tiny;
- return tiny-one;
+ return 1.0 - tiny;
+ return tiny - 1.0;
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
SET_LOW_WORD(z,0);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx >= 0)
sb5+s*(sb6+s*sb7))))));
}
z = x;
SET_LOW_WORD(z,0);
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx >= 0)
- return one-r/x;
- return r/x-one;
+ return 1.0 - r/x;
+ return r/x - 1.0;
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
- return (double)(((uint32_t)hx>>31)<<1) + one/x;
+ return (double)(((uint32_t)hx>>31)<<1) + 1.0/x;
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */
}
if (ix < 0x3feb0000) { /* |x| < 0.84375 */
if (ix < 0x3c700000) /* |x| < 2**-56 */
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+ s = 1.0+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if (hx < 0x3fd00000) { /* x < 1/4 */
y = r/s;
if (hx < 0x3fd00000) { /* x < 1/4 */
- r += x-half;
- return half - r ;
+ r += x - 0.5;
+ return 0.5 - r ;
}
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
}
}
if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+ Q = 1.0+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
return z - P/Q;
} else {
z = erx+P/Q;
return z - P/Q;
} else {
z = erx+P/Q;
}
}
if (ix < 0x403c0000) { /* |x| < 28 */
x = fabs(x);
}
}
if (ix < 0x403c0000) { /* |x| < 28 */
x = fabs(x);
if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = 1.0+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40180000) /* x < -6 */
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40180000) /* x < -6 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = 1.0+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
sb5+s*(sb6+s*sb7))))));
}
z = x;
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx > 0)
return r/x;
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if (hx > 0)
return r/x;
}
if (hx > 0)
return tiny*tiny;
}
if (hx > 0)
return tiny*tiny;
static const float
tiny = 1e-30,
static const float
tiny = 1e-30,
-half = 5.0000000000e-01, /* 0x3F000000 */
-one = 1.0000000000e+00, /* 0x3F800000 */
-two = 2.0000000000e+00, /* 0x40000000 */
/* c = (subfloat)0.84506291151 */
erx = 8.4506291151e-01, /* 0x3f58560b */
/*
/* c = (subfloat)0.84506291151 */
erx = 8.4506291151e-01, /* 0x3f58560b */
/*
if (ix >= 0x7f800000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
if (ix >= 0x7f800000) {
/* erf(nan)=nan, erf(+-inf)=+-1 */
i = ((uint32_t)hx>>31)<<1;
- return (float)(1-i)+one/x;
+ return (float)(1-i)+1.0f/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x31800000) { /* |x| < 2**-28 */
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x31800000) { /* |x| < 2**-28 */
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+ s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
y = r/s;
return x + x*y;
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+ Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40c00000) { /* inf > |x| >= 6 */
if (hx >= 0)
if (hx >= 0)
return erx + P/Q;
return -erx - P/Q;
}
if (ix >= 0x40c00000) { /* inf > |x| >= 6 */
if (hx >= 0)
- return one - tiny;
- return tiny - one;
+ return 1.0f - tiny;
+ return tiny - 1.0f;
if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
if (ix < 0x4036DB6E) { /* |x| < 1/0.35 */
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(z, ix&0xfffff000);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx >= 0)
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
SET_FLOAT_WORD(z, ix&0xfffff000);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx >= 0)
- return one - r/x;
- return r/x - one;
+ return 1.0f - r/x;
+ return r/x - 1.0f;
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000) {
/* erfc(nan)=nan, erfc(+-inf)=0,2 */
- return (float)(((uint32_t)hx>>31)<<1) + one/x;
+ return (float)(((uint32_t)hx>>31)<<1) + 1.0f/x;
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x23800000) /* |x| < 2**-56 */
}
if (ix < 0x3f580000) { /* |x| < 0.84375 */
if (ix < 0x23800000) /* |x| < 2**-56 */
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
- s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+ s = 1.0f+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if (hx < 0x3e800000) { /* x<1/4 */
y = r/s;
if (hx < 0x3e800000) { /* x<1/4 */
- r += (x-half);
- return half - r ;
+ r += (x-0.5f);
+ return 0.5f - r ;
}
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
}
}
if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
- Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
+ Q = 1.0f+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
return z - P/Q;
} else {
z = erx + P/Q;
return z - P/Q;
} else {
z = erx + P/Q;
}
}
if (ix < 0x41e00000) { /* |x| < 28 */
x = fabsf(x);
}
}
if (ix < 0x41e00000) { /* |x| < 28 */
x = fabsf(x);
if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
if (ix < 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = 1.0f+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40c00000) /* x < -6 */
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if (hx < 0 && ix >= 0x40c00000) /* x < -6 */
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = 1.0f+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
sb5+s*(sb6+s*sb7))))));
}
GET_FLOAT_WORD(ix, x);
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx > 0)
return r/x;
r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
if (hx > 0)
return r/x;
}
if (hx > 0)
return tiny*tiny;
}
if (hx > 0)
return tiny*tiny;
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
tiny = 1e-4931L,
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
tiny = 1e-4931L,
-half = 0.5L,
-one = 1.0L,
-two = 2.0L,
/* c = (float)0.84506291151 */
erx = 0.845062911510467529296875L,
/* c = (float)0.84506291151 */
erx = 0.845062911510467529296875L,
int32_t ix, i;
uint32_t se, i0, i1;
int32_t ix, i;
uint32_t se, i0, i1;
- GET_LDOUBLE_WORDS (se, i0, i1, x);
+ GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erf(nan)=nan */
i = ((se & 0xffff) >> 15) << 1;
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erf(nan)=nan */
i = ((se & 0xffff) >> 15) << 1;
- return (long double)(1 - i) + one / x; /* erf(+-inf)=+-1 */
+ return (long double)(1 - i) + 1.0 / x; /* erf(+-inf)=+-1 */
}
ix = (ix << 16) | (i0 >> 16);
}
ix = (ix << 16) | (i0 >> 16);
return x + x * y;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
return x + x * y;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
}
if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */
if ((se & 0x8000) == 0)
}
if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */
if ((se & 0x8000) == 0)
- return one - tiny;
- return tiny - one;
+ return 1.0 - tiny;
+ return tiny - 1.0;
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
SET_LDOUBLE_WORDS(z, i, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
SET_LDOUBLE_WORDS(z, i, i0, i1);
r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
- return one - r / x;
- return r / x - one;
+ return 1.0 - r / x;
+ return r / x - 1.0;
}
long double erfcl(long double x)
}
long double erfcl(long double x)
long double R, S, P, Q, s, y, z, r;
uint32_t se, i0, i1;
long double R, S, P, Q, s, y, z, r;
uint32_t se, i0, i1;
- GET_LDOUBLE_WORDS (se, i0, i1, x);
+ GET_LDOUBLE_WORDS(se, i0, i1, x);
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
ix = se & 0x7fff;
if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
- return (long double)(((se & 0xffff) >> 15) << 1) + one / x;
+ return (long double)(((se & 0xffff) >> 15) << 1) + 1.0 / x;
}
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fbe0000) /* |x| < 2**-65 */
}
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffed800) { /* |x| < 0.84375 */
if (ix < 0x3fbe0000) /* |x| < 2**-65 */
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
z = x * x;
r = pp[0] + z * (pp[1] +
z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
if (ix < 0x3ffd8000) /* x < 1/4 */
z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
y = r / s;
if (ix < 0x3ffd8000) /* x < 1/4 */
- return one - (x + x * y);
+ return 1.0 - (x + x * y);
- r += x - half;
- return half - r;
+ r += x - 0.5L;
+ return 0.5L - r;
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
}
if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0) {
P = pa[0] + s * (pa[1] + s * (pa[2] +
s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
Q = qa[0] + s * (qa[1] + s * (qa[2] +
s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
if ((se & 0x8000) == 0) {
return z - P / Q;
}
z = erx + P / Q;
return z - P / Q;
}
z = erx + P / Q;
}
if (ix < 0x4005d600) { /* |x| < 107 */
}
if (ix < 0x4005d600) { /* |x| < 107 */
- x = fabsl (x);
- s = one / (x * x);
+ x = fabsl(x);
+ s = 1.0 / (x * x);
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
s * (sb[5] + s * (sb[6] + s))))));
} else { /* 107 > |x| >= 6.666 */
if (se & 0x8000)
s * (sb[5] + s * (sb[6] + s))))));
} else { /* 107 > |x| >= 6.666 */
if (se & 0x8000)
- return two - tiny;/* x < -6.666 */
+ return 2.0 - tiny;/* x < -6.666 */
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
z = x;
R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
s * (rc[4] + s * rc[5]))));
S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
s * (sc[4] + s))));
}
z = x;
- GET_LDOUBLE_WORDS (hx, i0, i1, z);
+ GET_LDOUBLE_WORDS(hx, i0, i1, z);
i1 = 0;
i0 &= 0xffffff00;
i1 = 0;
i0 &= 0xffffff00;
- SET_LDOUBLE_WORDS (z, hx, i0, i1);
- r = expl (-z * z - 0.5625) *
- expl ((z - x) * (z + x) + R / S);
+ SET_LDOUBLE_WORDS(z, hx, i0, i1);
+ r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
if ((se & 0x8000) == 0)
return r / x;
if ((se & 0x8000) == 0)
return r / x;
}
if ((se & 0x8000) == 0)
return tiny * tiny;
}
if ((se & 0x8000) == 0)
return tiny * tiny;
#include "libm.h"
static const double
#include "libm.h"
static const double
halF[2] = {0.5,-0.5,},
huge = 1.0e+300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
halF[2] = {0.5,-0.5,},
huge = 1.0e+300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
STRICT_ASSIGN(double, x, hi - lo);
} else if(hx < 0x3e300000) { /* |x| < 2**-28 */
/* raise inexact */
STRICT_ASSIGN(double, x, hi - lo);
} else if(hx < 0x3e300000) { /* |x| < 2**-28 */
/* raise inexact */
- if (huge+x > one)
- return one+x;
+ if (huge+x > 1.0)
+ return 1.0+x;
INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0);
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if (k == 0)
INSERT_WORDS(twopk, 0x3ff00000+((k+1000)<<20), 0);
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if (k == 0)
- return one - ((x*c)/(c-2.0) - x);
- y = one-((lo-(x*c)/(2.0-c))-hi);
+ return 1.0 - ((x*c)/(c-2.0) - x);
+ y = 1.0-((lo-(x*c)/(2.0-c))-hi);
if (k < -1021)
return y*twopk*twom1000;
if (k == 1024)
if (k < -1021)
return y*twopk*twom1000;
if (k == 1024)
#include "libm.h"
static const float
#include "libm.h"
static const float
halF[2] = {0.5,-0.5,},
huge = 1.0e+30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
halF[2] = {0.5,-0.5,},
huge = 1.0e+30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
STRICT_ASSIGN(float, x, hi - lo);
} else if(hx < 0x39000000) { /* |x|<2**-14 */
/* raise inexact */
STRICT_ASSIGN(float, x, hi - lo);
} else if(hx < 0x39000000) { /* |x|<2**-14 */
/* raise inexact */
- if (huge+x > one)
- return one + x;
+ if (huge+x > 1.0f)
+ return 1.0f + x;
SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23));
c = x - t*(P1+t*P2);
if (k == 0)
SET_FLOAT_WORD(twopk, 0x3f800000+((k+100)<<23));
c = x - t*(P1+t*P2);
if (k == 0)
- return one - ((x*c)/(c - 2.0f) - x);
- y = one - ((lo - (x*c)/(2.0f - c)) - hi);
+ return 1.0f - ((x*c)/(c - 2.0f) - x);
+ y = 1.0f - ((lo - (x*c)/(2.0f - c)) - hi);
if (k < -125)
return y*twopk*twom100;
if (k == 128)
if (k < -125)
return y*twopk*twom100;
if (k == 128)
#include "libm.h"
static const double
#include "libm.h"
static const double
huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */
/* raise inexact */
if(x+tiny<0.0)
if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */
/* raise inexact */
if(x+tiny<0.0)
- return tiny-one; /* return -1 */
+ return tiny-1.0; /* return -1 */
/* x is now in primary range */
hfx = 0.5*x;
hxs = x*hfx;
/* x is now in primary range */
hfx = 0.5*x;
hxs = x*hfx;
- r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
+ r1 = 1.0+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
t = 3.0-r1*hfx;
e = hxs*((r1-t)/(6.0 - x*t));
if (k == 0) /* c is 0 */
t = 3.0-r1*hfx;
e = hxs*((r1-t)/(6.0 - x*t));
if (k == 0) /* c is 0 */
if (k == 1) {
if (x < -0.25)
return -2.0*(e-(x+0.5));
if (k == 1) {
if (x < -0.25)
return -2.0*(e-(x+0.5));
}
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
}
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
if (k == 1024)
y = y*2.0*0x1p1023;
else
y = y*twopk;
if (k == 1024)
y = y*2.0*0x1p1023;
else
y = y*twopk;
if (k < 20) {
SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
y = t-(e-x);
if (k < 20) {
SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
y = t-(e-x);
} else {
SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */
y = x-(e+t);
} else {
SET_HIGH_WORD(t, ((0x3ff-k)<<20)); /* 2^-k */
y = x-(e+t);
#include "libm.h"
static const float
#include "libm.h"
static const float
huge = 1.0e+30,
tiny = 1.0e-30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
huge = 1.0e+30,
tiny = 1.0e-30,
o_threshold = 8.8721679688e+01, /* 0x42b17180 */
if (xsb != 0) { /* x < -27*ln2 */
/* raise inexact */
if (x+tiny < 0.0f)
if (xsb != 0) { /* x < -27*ln2 */
/* raise inexact */
if (x+tiny < 0.0f)
- return tiny-one; /* return -1 */
+ return tiny-1.0f; /* return -1 */
/* x is now in primary range */
hfx = 0.5f*x;
hxs = x*hfx;
/* x is now in primary range */
hfx = 0.5f*x;
hxs = x*hfx;
- r1 = one+hxs*(Q1+hxs*Q2);
+ r1 = 1.0f+hxs*(Q1+hxs*Q2);
t = 3.0f - r1*hfx;
e = hxs*((r1-t)/(6.0f - x*t));
if (k == 0) /* c is 0 */
t = 3.0f - r1*hfx;
e = hxs*((r1-t)/(6.0f - x*t));
if (k == 0) /* c is 0 */
if (k == 1) {
if (x < -0.25f)
return -2.0f*(e-(x+0.5f));
if (k == 1) {
if (x < -0.25f)
return -2.0f*(e-(x+0.5f));
- return one + 2.0f*(x-e);
+ return 1.0f + 2.0f*(x-e);
}
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
}
if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */
if (k == 128)
y = y*2.0f*0x1p127f;
else
y = y*twopk;
if (k == 128)
y = y*2.0f*0x1p127f;
else
y = y*twopk;
if (k < 23) {
SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
y = t - (e - x);
if (k < 23) {
SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
y = t - (e - x);
} else {
SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */
y = x - (e + t);
} else {
SET_FLOAT_WORD(t, (0x7f-k)<<23); /* 2^-k */
y = x - (e + t);
static const double
huge = 1e300,
static const double
huge = 1e300,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0, 2.00] */
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0, 2.00] */
S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
S03 = 5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
-static const double zero = 0.0;
-
double j0(double x)
{
double z, s,c,ss,cc,r,u,v;
double j0(double x)
{
double z, s,c,ss,cc,r,u,v;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
x = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
cc = s+c;
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
cc = s+c;
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
}
if (ix < 0x3f200000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
}
if (ix < 0x3f200000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
if (ix < 0x3e400000) /* |x| < 2**-27 */
if (ix < 0x3e400000) /* |x| < 2**-27 */
- return one;
- return one - 0.25*x*x;
+ return 1.0;
+ return 1.0 - 0.25*x*x;
}
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
}
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
- s = one+z*(S01+z*(S02+z*(S03+z*S04)));
+ s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
if (ix < 0x3FF00000) { /* |x| < 1.00 */
if (ix < 0x3FF00000) { /* |x| < 1.00 */
- return one + z*(-0.25+(r/s));
+ return 1.0 + z*(-0.25+(r/s));
- return (one+u)*(one-u) + z*(r/s);
+ return (1.0+u)*(1.0-u) + z*(r/s);
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7ff00000)
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7ff00000)
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
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
*/
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
*/
if (ix < 0x7fe00000) { /* make sure x+x does not overflow */
z = -cos(x+x);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
- v = one+z*(v01+z*(v02+z*(v03+z*v04)));
+ v = 1.0+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0(x)*log(x));
}
return u/v + tpi*(j0(x)*log(x));
}
else if (ix >= 0x40122E8B){p = pR5; q = pS5;}
else if (ix >= 0x4006DB6D){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
else if (ix >= 0x40122E8B){p = pR5; q = pS5;}
else if (ix >= 0x4006DB6D){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
- return one + r/s;
+ s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+ return 1.0 + r/s;
else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+ s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125 + r/s)/x;
}
return (-.125 + r/s)/x;
}
static const float
huge = 1e30,
static const float
huge = 1e30,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0, 2.00] */
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0, 2.00] */
S03 = 5.1354652442e-07, /* 0x3509daa6 */
S04 = 1.1661400734e-09; /* 0x30a045e8 */
S03 = 5.1354652442e-07, /* 0x3509daa6 */
S04 = 1.1661400734e-09; /* 0x30a045e8 */
-static const float zero = 0.0;
-
float j0f(float x)
{
float z, s,c,ss,cc,r,u,v;
float j0f(float x)
{
float z, s,c,ss,cc,r,u,v;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000)
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000)
x = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
x = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
cc = s+c;
if (ix < 0x7f000000) { /* make sure x+x does not overflow */
z = -cosf(x+x);
cc = s+c;
if (ix < 0x7f000000) { /* make sure x+x does not overflow */
z = -cosf(x+x);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
}
if (ix < 0x39000000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
}
if (ix < 0x39000000) { /* |x| < 2**-13 */
/* raise inexact if x != 0 */
if (ix < 0x32000000) /* |x| < 2**-27 */
if (ix < 0x32000000) /* |x| < 2**-27 */
- return one;
- return one - 0.25f*x*x;
+ return 1.0f;
+ return 1.0f - 0.25f*x*x;
}
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
}
}
z = x*x;
r = z*(R02+z*(R03+z*(R04+z*R05)));
- s = one+z*(S01+z*(S02+z*(S03+z*S04)));
+ s = 1.0f+z*(S01+z*(S02+z*(S03+z*S04)));
if(ix < 0x3F800000) { /* |x| < 1.00 */
if(ix < 0x3F800000) { /* |x| < 1.00 */
- return one + z*(-0.25f + (r/s));
+ return 1.0f + z*(-0.25f + (r/s));
- return (one+u)*(one-u) + z*(r/s);
+ return (1.0f+u)*(1.0f-u) + z*(r/s);
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7f800000)
ix = 0x7fffffff & hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if (ix >= 0x7f800000)
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
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
*/
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = -cosf(x+x);
*/
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = -cosf(x+x);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
}
z = x*x;
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
- v = one+z*(v01+z*(v02+z*(v03+z*v04)));
+ v = 1.0f+z*(v01+z*(v02+z*(v03+z*v04)));
return u/v + tpi*(j0f(x)*logf(x));
}
return u/v + tpi*(j0f(x)*logf(x));
}
else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
else if (ix >= 0x4036db68){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
else if (ix >= 0x40f71c58){p = pR5; q = pS5;}
else if (ix >= 0x4036db68){p = pR3; q = pS3;}
else if (ix >= 0x40000000){p = pR2; q = pS2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
- return one + r/s;
+ s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+ return 1.0f + r/s;
else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
else if (ix >= 0x4036db68){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
else if (ix >= 0x4036db68){p = qR3; q = qS3;}
else if (ix >= 0x40000000){p = qR2; q = qS2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+ s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (-.125f + r/s)/x;
}
return (-.125f + r/s)/x;
}
static const double
huge = 1e300,
static const double
huge = 1e300,
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0,2] */
invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
/* R0/S0 on [0,2] */
s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
-static const double zero = 0.0;
-
double j1(double x)
{
double z,s,c,ss,cc,r,u,v,y;
double j1(double x)
{
double z,s,c,ss,cc,r,u,v,y;
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
GET_HIGH_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000)
y = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(y);
y = fabs(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(y);
cc = s-c;
if (ix < 0x7fe00000) { /* make sure y+y not overflow */
z = cos(y+y);
cc = s-c;
if (ix < 0x7fe00000) { /* make sure y+y not overflow */
z = cos(y+y);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
}
if (ix < 0x3e400000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */
}
if (ix < 0x3e400000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */
return 0.5*x;
}
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
return 0.5*x;
}
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
- s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
+ s = 1.0+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
r *= x;
return x*0.5 + r/s;
}
r *= x;
return x*0.5 + r/s;
}
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7ff00000)
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7ff00000)
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
cc = s-c;
if (ix < 0x7fe00000) { /* make sure x+x not overflow */
z = cos(x+x);
cc = s-c;
if (ix < 0x7fe00000) { /* make sure x+x not overflow */
z = cos(x+x);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
- v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
- return x*(u/v) + tpi*(j1(x)*log(x)-one/x);
+ v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
+ return x*(u/v) + tpi*(j1(x)*log(x)-1.0/x);
}
/* For x >= 8, the asymptotic expansions of pone is
}
/* For x >= 8, the asymptotic expansions of pone is
else if (ix >= 0x40122E8B){p = pr5; q = ps5;}
else if (ix >= 0x4006DB6D){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;}
else if (ix >= 0x40122E8B){p = pr5; q = ps5;}
else if (ix >= 0x4006DB6D){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
- return one+ r/s;
+ s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+ return 1.0+ r/s;
}
/* For x >= 8, the asymptotic expansions of qone is
}
/* For x >= 8, the asymptotic expansions of qone is
else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;}
else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+ s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
static const float
huge = 1e30,
static const float
huge = 1e30,
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0,2] */
invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
tpi = 6.3661974669e-01, /* 0x3f22f983 */
/* R0/S0 on [0,2] */
s04 = 5.0463624390e-09, /* 0x31ad6446 */
s05 = 1.2354227016e-11; /* 0x2d59567e */
s04 = 5.0463624390e-09, /* 0x31ad6446 */
s05 = 1.2354227016e-11; /* 0x2d59567e */
-static const float zero = 0.0;
-
float j1f(float x)
{
float z,s,c,ss,cc,r,u,v,y;
float j1f(float x)
{
float z,s,c,ss,cc,r,u,v,y;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000)
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix >= 0x7f800000)
y = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(y);
y = fabsf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(y);
cc = s-c;
if (ix < 0x7f000000) { /* make sure y+y not overflow */
z = cosf(y+y);
cc = s-c;
if (ix < 0x7f000000) { /* make sure y+y not overflow */
z = cosf(y+y);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
}
if (ix < 0x32000000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */
}
if (ix < 0x32000000) { /* |x| < 2**-27 */
/* raise inexact if x!=0 */
return 0.5f*x;
}
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
return 0.5f*x;
}
z = x*x;
r = z*(r00+z*(r01+z*(r02+z*r03)));
- s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
+ s = 1.0f+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
r *= x;
return 0.5f*x + r/s;
}
r *= x;
return 0.5f*x + r/s;
}
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7f800000)
ix = 0x7fffffff & hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if (ix >= 0x7f800000)
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
if (ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
cc = s-c;
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = cosf(x+x);
cc = s-c;
if (ix < 0x7f000000) { /* make sure x+x not overflow */
z = cosf(x+x);
cc = z/ss;
else
ss = z/cc;
cc = z/ss;
else
ss = z/cc;
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
return -tpi/x;
z = x*x;
u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
- v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
- return x*(u/v) + tpi*(j1f(x)*logf(x)-one/x);
+ v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
+ return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x);
}
/* For x >= 8, the asymptotic expansions of pone is
}
/* For x >= 8, the asymptotic expansions of pone is
else if (ix >= 0x40f71c58){p = pr5; q = ps5;}
else if (ix >= 0x4036db68){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;}
else if (ix >= 0x40f71c58){p = pr5; q = ps5;}
else if (ix >= 0x4036db68){p = pr3; q = ps3;}
else if (ix >= 0x40000000){p = pr2; q = ps2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
- return one + r/s;
+ s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+ return 1.0f + r/s;
}
/* For x >= 8, the asymptotic expansions of qone is
}
/* For x >= 8, the asymptotic expansions of qone is
else if (ix >= 0x40f71c58){p = qr5; q = qs5;}
else if (ix >= 0x4036db68){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;}
else if (ix >= 0x40f71c58){p = qr5; q = qs5;}
else if (ix >= 0x4036db68){p = qr3; q = qs3;}
else if (ix >= 0x40000000){p = qr2; q = qs2;}
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
- s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+ s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
return (.375f + r/s)/x;
}
return (.375f + r/s)/x;
}
-static const double
-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 invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
double jn(int n, double x)
{
double jn(int n, double x)
{
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabs(x);
if ((ix|lx) == 0 || ix >= 0x7ff00000) /* if x is 0 or inf */
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabs(x);
if ((ix|lx) == 0 || ix >= 0x7ff00000) /* if x is 0 or inf */
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 */
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 */
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
else {
temp = x*0.5;
b = temp;
else {
temp = x*0.5;
b = temp;
- for (a=one,i=2; i<=n; i++) {
+ for (a=1.0,i=2; i<=n; i++) {
a *= (double)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
a *= (double)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
- for (t=zero, i = 2*(n+k); i>=m; i -= 2)
- t = one/(i/x-t);
+ for (t=0.0, i = 2*(n+k); i>=m; i -= 2)
+ t = 1.0/(i/x-t);
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* likely underflow to zero
*/
tmp = n;
* likely underflow to zero
*/
tmp = n;
tmp = tmp*log(fabs(v*tmp));
if (tmp < 7.09782712893383973096e+02) {
for (i=n-1,di=(double)(i+i); i>0; i--) {
tmp = tmp*log(fabs(v*tmp));
if (tmp < 7.09782712893383973096e+02) {
for (i=n-1,di=(double)(i+i); i>0; i--) {
b *= di;
b = b/x - a;
a = temp;
b *= di;
b = b/x - a;
a = temp;
}
} else {
for (i=n-1,di=(double)(i+i); i>0; i--) {
}
} else {
for (i=n-1,di=(double)(i+i); i>0; i--) {
b *= di;
b = b/x - a;
a = temp;
b *= di;
b = b/x - a;
a = temp;
/* scale b to avoid spurious overflow */
if (b > 1e100) {
a /= b;
t /= b;
/* scale b to avoid spurious overflow */
if (b > 1e100) {
a /= b;
t /= b;
if ((ix|((uint32_t)(lx|-lx))>>31) > 0x7ff00000)
return x+x;
if ((ix|lx) == 0)
if ((ix|((uint32_t)(lx|-lx))>>31) > 0x7ff00000)
return x+x;
if ((ix|lx) == 0)
sign = 1;
if (n < 0) {
n = -n;
sign = 1;
if (n < 0) {
n = -n;
if (n == 1)
return sign*y1(x);
if (ix == 0x7ff00000)
if (n == 1)
return sign*y1(x);
if (ix == 0x7ff00000)
if (ix >= 0x52D00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
if (ix >= 0x52D00000) { /* x > 2**302 */
/* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
#define _GNU_SOURCE
#include "libm.h"
#define _GNU_SOURCE
#include "libm.h"
-static const float
-two = 2.0000000000e+00, /* 0x40000000 */
-one = 1.0000000000e+00; /* 0x3F800000 */
-
-static const float zero = 0.0000000000e+00;
-
float jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
float jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
else if((float)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x);
else if((float)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x);
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
else {
temp = 0.5f * x;
b = temp;
else {
temp = 0.5f * x;
b = temp;
- for (a=one,i=2; i<=n; i++) {
+ for (a=1.0f,i=2; i<=n; i++) {
a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
- for (t=zero, i = 2*(n+k); i>=m; i -= 2)
- t = one/(i/x-t);
+ for (t=0.0f, i = 2*(n+k); i>=m; i -= 2)
+ t = 1.0f/(i/x-t);
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* likely underflow to zero
*/
tmp = n;
* likely underflow to zero
*/
tmp = n;
tmp = tmp*logf(fabsf(v*tmp));
if (tmp < 88.721679688f) {
for (i=n-1,di=(float)(i+i); i>0; i--) {
tmp = tmp*logf(fabsf(v*tmp));
if (tmp < 88.721679688f) {
for (i=n-1,di=(float)(i+i); i>0; i--) {
b *= di;
b = b/x - a;
a = temp;
b *= di;
b = b/x - a;
a = temp;
}
} else {
for (i=n-1,di=(float)(i+i); i>0; i--){
}
} else {
for (i=n-1,di=(float)(i+i); i>0; i--){
b *= di;
b = b/x - a;
a = temp;
b *= di;
b = b/x - a;
a = temp;
/* scale b to avoid spurious overflow */
if (b > 1e10f) {
a /= b;
t /= b;
/* scale b to avoid spurious overflow */
if (b > 1e10f) {
a /= b;
t /= b;
if (ix > 0x7f800000)
return x+x;
if (ix == 0)
if (ix > 0x7f800000)
return x+x;
if (ix == 0)
sign = 1;
if (n < 0) {
n = -n;
sign = 1;
if (n < 0) {
n = -n;
if (n == 1)
return sign*y1f(x);
if (ix == 0x7f800000)
if (n == 1)
return sign*y1f(x);
if (ix == 0x7f800000)
static const double
two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
static const double
two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */
a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
-static const double zero = 0.00000000000000000000e+00;
-
static double sin_pi(double x)
{
double y,z;
static double sin_pi(double x)
{
double y,z;
ix &= 0x7fffffff;
if (ix < 0x3fd00000)
ix &= 0x7fffffff;
if (ix < 0x3fd00000)
- return __sin(pi*x, zero, 0);
+ return __sin(pi*x, 0.0, 0);
y = -x; /* negative x is assumed */
y = -x; /* negative x is assumed */
n = (int)(y*4.0);
} else {
if (ix >= 0x43400000) {
n = (int)(y*4.0);
} else {
if (ix >= 0x43400000) {
- y = zero; /* y must be even */
+ y = 0.0; /* y must be even */
n = 0;
} else {
if (ix < 0x43300000)
n = 0;
} else {
if (ix < 0x43300000)
- case 0: y = __sin(pi*y, zero, 0); break;
+ case 0: y = __sin(pi*y, 0.0, 0); break;
- case 2: y = __cos(pi*(0.5-y), zero); break;
+ case 2: y = __cos(pi*(0.5-y), 0.0); break;
- case 4: y = __sin(pi*(one-y), zero, 0); break;
+ case 4: y = __sin(pi*(1.0-y), 0.0, 0); break;
- case 6: y = -__cos(pi*(y-1.5), zero); break;
- default: y = __sin(pi*(y-2.0), zero, 0); break;
+ case 6: y = -__cos(pi*(y-1.5), 0.0); break;
+ default: y = __sin(pi*(y-2.0), 0.0, 0); break;
if (ix >= 0x7ff00000)
return x*x;
if ((ix|lx) == 0)
if (ix >= 0x7ff00000)
return x*x;
if ((ix|lx) == 0)
if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx < 0) {
*signgamp = -1;
if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx < 0) {
*signgamp = -1;
}
if (hx < 0) {
if (ix >= 0x43300000) /* |x|>=2**52, must be -integer */
}
if (hx < 0) {
if (ix >= 0x43300000) /* |x|>=2**52, must be -integer */
- if (t == zero) /* -integer */
- return one/zero;
+ if (t == 0.0) /* -integer */
+ return 1.0/0.0;
nadj = log(pi/fabs(t*x));
nadj = log(pi/fabs(t*x));
*signgamp = -1;
x = -x;
}
*signgamp = -1;
x = -x;
}
if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -log(x);
if (ix >= 0x3FE76944) {
if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -log(x);
if (ix >= 0x3FE76944) {
i = 0;
} else if (ix >= 0x3FCDA661) {
i = 0;
} else if (ix >= 0x3FCDA661) {
i = 1;
} else {
y = x;
i = 2;
}
} else {
i = 1;
} else {
y = x;
i = 2;
}
} else {
if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */
y = 2.0 - x;
i = 0;
if (ix >= 0x3FFBB4C3) { /* [1.7316,2] */
y = 2.0 - x;
i = 0;
y = x - tc;
i = 1;
} else {
y = x - tc;
i = 1;
} else {
break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
- p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
+ p2 = 1.0+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5*y + p1/p2;
}
} else if (ix < 0x40200000) { /* x < 8.0 */
i = (int)x;
y = x - (double)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
r += -0.5*y + p1/p2;
}
} else if (ix < 0x40200000) { /* x < 8.0 */
i = (int)x;
y = x - (double)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
- q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
- r = half*y+p/q;
- z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
+ q = 1.0+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
+ r = 0.5*y+p/q;
+ z = 1.0; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7: z *= y + 6.0; /* FALLTHRU */
case 6: z *= y + 5.0; /* FALLTHRU */
switch (i) {
case 7: z *= y + 6.0; /* FALLTHRU */
case 6: z *= y + 5.0; /* FALLTHRU */
}
} else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */
t = log(x);
}
} else if (ix < 0x43900000) { /* 8.0 <= x < 2**58 */
t = log(x);
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
- r = (x-half)*(t-one)+w;
} else /* 2**58 <= x <= inf */
} else /* 2**58 <= x <= inf */
if (hx < 0)
r = nadj - r;
return r;
if (hx < 0)
r = nadj - r;
return r;
static const float
two23= 8.3886080000e+06, /* 0x4b000000 */
static const float
two23= 8.3886080000e+06, /* 0x4b000000 */
-half= 5.0000000000e-01, /* 0x3f000000 */
-one = 1.0000000000e+00, /* 0x3f800000 */
pi = 3.1415927410e+00, /* 0x40490fdb */
a0 = 7.7215664089e-02, /* 0x3d9e233f */
a1 = 3.2246702909e-01, /* 0x3ea51a66 */
pi = 3.1415927410e+00, /* 0x40490fdb */
a0 = 7.7215664089e-02, /* 0x3d9e233f */
a1 = 3.2246702909e-01, /* 0x3ea51a66 */
w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
-static const float zero = 0.0000000000e+00;
-
static float sin_pif(float x)
{
float y,z;
static float sin_pif(float x)
{
float y,z;
n = (int)(y*4.0f);
} else {
if (ix >= 0x4b800000) {
n = (int)(y*4.0f);
} else {
if (ix >= 0x4b800000) {
- y = zero; /* y must be even */
+ y = 0.0f; /* y must be even */
n = 0;
} else {
if (ix < 0x4b000000)
n = 0;
} else {
if (ix < 0x4b000000)
case 1:
case 2: y = __cosdf(pi*(0.5f - y)); break;
case 3:
case 1:
case 2: y = __cosdf(pi*(0.5f - y)); break;
case 3:
- case 4: y = __sindf(pi*(one - y)); break;
+ case 4: y = __sindf(pi*(1.0f - y)); break;
case 5:
case 6: y = -__cosdf(pi*(y - 1.5f)); break;
default: y = __sindf(pi*(y - 2.0f)); break;
case 5:
case 6: y = -__cosdf(pi*(y - 1.5f)); break;
default: y = __sindf(pi*(y - 2.0f)); break;
if (ix >= 0x7f800000)
return x*x;
if (ix == 0)
if (ix >= 0x7f800000)
return x*x;
if (ix == 0)
if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */
if (hx < 0) {
*signgamp = -1;
if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */
if (hx < 0) {
*signgamp = -1;
}
if (hx < 0) {
if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */
}
if (hx < 0) {
if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */
- if (t == zero) /* -integer */
- return one/zero;
+ if (t == 0.0f) /* -integer */
+ return 1.0f/0.0f;
nadj = logf(pi/fabsf(t*x));
nadj = logf(pi/fabsf(t*x));
*signgamp = -1;
x = -x;
}
*signgamp = -1;
x = -x;
}
if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -logf(x);
if (ix >= 0x3f3b4a20) {
if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
r = -logf(x);
if (ix >= 0x3f3b4a20) {
i = 0;
} else if (ix >= 0x3e6d3308) {
i = 0;
} else if (ix >= 0x3e6d3308) {
i = 1;
} else {
y = x;
i = 2;
}
} else {
i = 1;
} else {
y = x;
i = 2;
}
} else {
if (ix >= 0x3fdda618) { /* [1.7316,2] */
y = 2.0f - x;
i = 0;
if (ix >= 0x3fdda618) { /* [1.7316,2] */
y = 2.0f - x;
i = 0;
y = x - tc;
i = 1;
} else {
y = x - tc;
i = 1;
} else {
break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
break;
case 2:
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
- p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
+ p2 = 1.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
r += -0.5f*y + p1/p2;
}
} else if (ix < 0x41000000) { /* x < 8.0 */
i = (int)x;
y = x - (float)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
r += -0.5f*y + p1/p2;
}
} else if (ix < 0x41000000) { /* x < 8.0 */
i = (int)x;
y = x - (float)i;
p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
- q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
- r = half*y+p/q;
- z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
+ q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
+ r = 0.5f*y+p/q;
+ z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7: z *= y + 6.0f; /* FALLTHRU */
case 6: z *= y + 5.0f; /* FALLTHRU */
switch (i) {
case 7: z *= y + 6.0f; /* FALLTHRU */
case 6: z *= y + 5.0f; /* FALLTHRU */
}
} else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */
t = logf(x);
}
} else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */
t = logf(x);
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
- r = (x-half)*(t-one)+w;
+ r = (x-0.5f)*(t-1.0f)+w;
} else /* 2**58 <= x <= inf */
} else /* 2**58 <= x <= inf */
if (hx < 0)
r = nadj - r;
return r;
if (hx < 0)
r = nadj - r;
return r;
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
-half = 0.5L,
-one = 1.0L,
pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L,
pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L,
w6 = -1.880801938119376907179E-3L,
w7 = 4.885026142432270781165E-3L;
w6 = -1.880801938119376907179E-3L,
w7 = 4.885026142432270781165E-3L;
-static const long double zero = 0.0L;
-
static long double sin_pi(long double x)
{
long double y, z;
static long double sin_pi(long double x)
{
long double y, z;
n = (int) (y*4.0);
} else {
if (ix >= 0x403f8000) { /* 2^64 */
n = (int) (y*4.0);
} else {
if (ix >= 0x403f8000) { /* 2^64 */
- y = zero; /* y must be even */
+ y = 0.0; /* y must be even */
n = 0;
} else {
if (ix < 0x403e8000) /* 2^63 */
n = 0;
} else {
if (ix < 0x403e8000) /* 2^63 */
- y = cosl(pi * (half - y));
+ y = cosl(pi * (0.5 - y));
- y = sinl(pi * (one - y));
+ y = sinl(pi * (1.0 - y));
if ((ix | i0 | i1) == 0) {
if (se & 0x8000)
*sg = -1;
if ((ix | i0 | i1) == 0) {
if (se & 0x8000)
*sg = -1;
}
ix = (ix << 16) | (i0 >> 16);
}
ix = (ix << 16) | (i0 >> 16);
}
if (se & 0x8000) {
t = sin_pi (x);
}
if (se & 0x8000) {
t = sin_pi (x);
- if (t == zero)
- return one / fabsl(t); /* -integer */
+ if (t == 0.0)
+ return 1.0 / fabsl(t); /* -integer */
nadj = logl(pi / fabsl(t * x));
nadj = logl(pi / fabsl(t * x));
else if (ix < 0x40008000) { /* x < 2.0 */
if (ix <= 0x3ffee666) { /* 8.99993896484375e-1 */
/* lgamma(x) = lgamma(x+1) - log(x) */
else if (ix < 0x40008000) { /* x < 2.0 */
if (ix <= 0x3ffee666) { /* 8.99993896484375e-1 */
/* lgamma(x) = lgamma(x+1) - log(x) */
if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */
if (ix >= 0x3ffebb4a) { /* 7.31597900390625e-1 */
i = 0;
} else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */
i = 0;
} else if (ix >= 0x3ffced33) { /* 2.31639862060546875e-1 */
i = 1;
} else { /* x < 0.23 */
y = x;
i = 2;
}
} else {
i = 1;
} else { /* x < 0.23 */
y = x;
i = 2;
}
} else {
if (ix >= 0x3fffdda6) { /* 1.73162841796875 */
/* [1.7316,2] */
y = x - 2.0;
if (ix >= 0x3fffdda6) { /* 1.73162841796875 */
/* [1.7316,2] */
y = x - 2.0;
i = 1;
} else {
/* [0.9, 1.23] */
i = 1;
} else {
/* [0.9, 1.23] */
case 0:
p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
case 0:
p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
- r += half * y + y * p1/p2;
+ r += 0.5 * y + y * p1/p2;
break;
case 1:
p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
break;
case 1:
p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
case 2:
p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
case 2:
p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
- r += (-half * y + p1 / p2);
+ r += (-0.5 * y + p1 / p2);
}
} else if (ix < 0x40028000) { /* 8.0 */
/* x < 8.0 */
i = (int)x;
}
} else if (ix < 0x40028000) { /* 8.0 */
/* x < 8.0 */
i = (int)x;
y = x - (double)i;
p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
y = x - (double)i;
p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
- r = half * y + p / q;
- z = one;/* lgamma(1+s) = log(s) + lgamma(s) */
+ r = 0.5 * y + p / q;
+ z = 1.0;/* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
case 7:
z *= (y + 6.0); /* FALLTHRU */
switch (i) {
case 7:
z *= (y + 6.0); /* FALLTHRU */
z *= (y + 3.0); /* FALLTHRU */
case 3:
z *= (y + 2.0); /* FALLTHRU */
z *= (y + 3.0); /* FALLTHRU */
case 3:
z *= (y + 2.0); /* FALLTHRU */
break;
}
} else if (ix < 0x40418000) { /* 2^66 */
/* 8.0 <= x < 2**66 */
break;
}
} else if (ix < 0x40418000) { /* 2^66 */
/* 8.0 <= x < 2**66 */
- t = logl (x);
- z = one / x;
+ t = logl(x);
+ z = 1.0 / x;
y = z * z;
w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
y = z * z;
w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
- r = (x - half) * (t - one) + w;
+ r = (x - 0.5) * (t - 1.0) + w;
} else /* 2**66 <= x <= inf */
} else /* 2**66 <= x <= inf */
- r = x * (logl (x) - one);
+ r = x * (logl(x) - 1.0);
if (se & 0x8000)
r = nadj - r;
return r;
if (se & 0x8000)
r = nadj - r;
return r;
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
-static const double zero = 0.0;
-
double log(double x)
{
double hfsq,f,s,z,R,w,t1,t2,dk;
double log(double x)
{
double hfsq,f,s,z,R,w,t1,t2,dk;
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
- return -two54/zero; /* log(+-0)=-inf */
+ return -two54/0.0; /* log(+-0)=-inf */
- return (x-x)/zero; /* log(-#) = NaN */
+ return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
/* subnormal number, scale up x */
k -= 54;
x *= two54;
k += i>>20;
f = x - 1.0;
if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */
k += i>>20;
f = x - 1.0;
if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */
}
dk = (double)k;
return dk*ln2_hi + dk*ln2_lo;
}
dk = (double)k;
return dk*ln2_hi + dk*ln2_lo;
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
-static const double zero = 0.0;
-
double log10(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
double log10(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
- return -two54/zero; /* log(+-0)=-inf */
+ return -two54/0.0; /* log(+-0)=-inf */
- return (x-x)/zero; /* log(-#) = NaN */
+ return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
/* subnormal number, scale up x */
k -= 54;
x *= two54;
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
- return zero; /* log(1) = +0 */
+ return 0.0; /* log(1) = +0 */
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */
log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */
-static const float zero = 0.0;
-
float log10f(float x)
{
float f,hfsq,hi,lo,r,y;
float log10f(float x)
{
float f,hfsq,hi,lo,r,y;
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
- return -two25/zero; /* log(+-0)=-inf */
+ return -two25/0.0f; /* log(+-0)=-inf */
- return (x-x)/zero; /* log(-#) = NaN */
+ return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
/* subnormal number, scale up x */
k -= 25;
x *= two25;
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
- return zero; /* log(1) = +0 */
+ return 0.0f; /* log(1) = +0 */
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;
Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
-static const double zero = 0.0;
-
double log1p(double x)
{
double hfsq,f,c,s,z,R,u;
double log1p(double x)
{
double hfsq,f,c,s,z,R,u;
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3ff00000) { /* x <= -1.0 */
if (x == -1.0)
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3ff00000) { /* x <= -1.0 */
if (x == -1.0)
- return -two54/zero; /* log1p(-1)=+inf */
+ return -two54/0.0; /* log1p(-1)=+inf */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if (ax < 0x3e200000) { /* |x| < 2**-29 */
/* raise inexact */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if (ax < 0x3e200000) { /* |x| < 2**-29 */
/* raise inexact */
- if (two54 + x > zero && ax < 0x3c900000) /* |x| < 2**-54 */
+ if (two54 + x > 0.0 && ax < 0x3c900000) /* |x| < 2**-54 */
return x;
return x - x*x*0.5;
}
return x;
return x - x*x*0.5;
}
}
hfsq = 0.5*f*f;
if (hu == 0) { /* |f| < 2**-20 */
}
hfsq = 0.5*f*f;
if (hu == 0) { /* |f| < 2**-20 */
c += k*ln2_lo;
return k*ln2_hi + c;
}
c += k*ln2_lo;
return k*ln2_hi + c;
}
Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */
Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */
-static const float zero = 0.0;
-
float log1pf(float x)
{
float hfsq,f,c,s,z,R,u;
float log1pf(float x)
{
float hfsq,f,c,s,z,R,u;
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3f800000) { /* x <= -1.0 */
if (x == -1.0f)
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
if (ax >= 0x3f800000) { /* x <= -1.0 */
if (x == -1.0f)
- return -two25/zero; /* log1p(-1)=+inf */
+ return -two25/0.0f; /* log1p(-1)=+inf */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if (ax < 0x38000000) { /* |x| < 2**-15 */
/* raise inexact */
return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if (ax < 0x38000000) { /* |x| < 2**-15 */
/* raise inexact */
- if (two25 + x > zero && ax < 0x33800000) /* |x| < 2**-24 */
+ if (two25 + x > 0.0f && ax < 0x33800000) /* |x| < 2**-24 */
return x;
return x - x*x*0.5f;
}
return x;
return x - x*x*0.5f;
}
}
hfsq = 0.5f * f * f;
if (hu == 0) { /* |f| < 2**-20 */
}
hfsq = 0.5f * f * f;
if (hu == 0) { /* |f| < 2**-20 */
c += k*ln2_lo;
return k*ln2_hi+c;
}
c += k*ln2_lo;
return k*ln2_hi+c;
}
if (xm1 == 0.0)
return xm1;
if (xm1 == 0.0)
return xm1;
/* Test for domain errors. */
/* Test for domain errors. */
- if (x <= 0.0L) {
- if (x == 0.0L)
+ if (x <= 0.0) {
+ if (x == 0.0)
return -INFINITY;
return NAN;
}
return -INFINITY;
return NAN;
}
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
- z = x - 0.5L;
- y = 0.5L * z + 0.5L;
+ z = x - 0.5;
+ y = 0.5 * z + 0.5;
} else { /* 2 (x-1)/(x+1) */
} else { /* 2 (x-1)/(x+1) */
- z = x - 0.5L;
- z -= 0.5L;
- y = 0.5L * x + 0.5L;
+ z = x - 0.5;
+ z -= 0.5;
+ y = 0.5 * x + 0.5;
if (x < SQRTH) {
e -= 1;
if (e != 0)
if (x < SQRTH) {
e -= 1;
if (e != 0)
else
x = xm1;
} else {
if (e != 0)
else
x = xm1;
} else {
if (e != 0)
ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
-static const double zero = 0.0;
-
double log2(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
double log2(double x)
{
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
k = 0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx) == 0)
- return -two54/zero; /* log(+-0)=-inf */
+ return -two54/0.0; /* log(+-0)=-inf */
- return (x-x)/zero; /* log(-#) = NaN */
+ return (x-x)/0.0; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 54;
x *= two54;
/* subnormal number, scale up x */
k -= 54;
x *= two54;
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
if (hx >= 0x7ff00000)
return x+x;
if (hx == 0x3ff00000 && lx == 0)
- return zero; /* log(1) = +0 */
+ return 0.0; /* log(1) = +0 */
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64) & 0x100000;
k += (hx>>20) - 1023;
hx &= 0x000fffff;
i = (hx+0x95f64) & 0x100000;
ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
-static const float zero = 0.0;
-
float log2f(float x)
{
float f,hfsq,hi,lo,r,y;
float log2f(float x)
{
float f,hfsq,hi,lo,r,y;
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
k = 0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff) == 0)
- return -two25/zero; /* log(+-0)=-inf */
+ return -two25/0.0f; /* log(+-0)=-inf */
- return (x-x)/zero; /* log(-#) = NaN */
+ return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
/* subnormal number, scale up x */
k -= 25;
x *= two25;
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
if (hx >= 0x7f800000)
return x+x;
if (hx == 0x3f800000)
- return zero; /* log(1) = +0 */
+ return 0.0f; /* log(1) = +0 */
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;
k += (hx>>23) - 127;
hx &= 0x007fffff;
i = (hx+(0x4afb0d))&0x800000;
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
-static const float zero = 0.0;
-
float logf(float x)
{
float hfsq,f,s,z,R,w,t1,t2,dk;
float logf(float x)
{
float hfsq,f,s,z,R,w,t1,t2,dk;
k = 0;
if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix & 0x7fffffff) == 0)
k = 0;
if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix & 0x7fffffff) == 0)
- return -two25/zero; /* log(+-0)=-inf */
+ return -two25/0.0f; /* log(+-0)=-inf */
- return (x-x)/zero; /* log(-#) = NaN */
+ return (x-x)/0.0f; /* log(-#) = NaN */
/* subnormal number, scale up x */
k -= 25;
x *= two25;
/* subnormal number, scale up x */
k -= 25;
x *= two25;
k += i>>23;
f = x - 1.0f;
if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */
k += i>>23;
f = x - 1.0f;
if ((0x007fffff & (0x8000 + ix)) < 0xc000) { /* -2**-9 <= f < 2**-9 */
dk = (float)k;
return dk*ln2_hi + dk*ln2_lo;
}
dk = (float)k;
return dk*ln2_hi + dk*ln2_lo;
}
-static const double one = 1.0;
-
double modf(double x, double *iptr)
{
int32_t i0,i1,j0;
double modf(double x, double *iptr)
{
int32_t i0,i1,j0;
*iptr = x;
return 0.0 / x;
}
*iptr = x;
return 0.0 / x;
}
GET_HIGH_WORD(high, x);
INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
return x;
GET_HIGH_WORD(high, x);
INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
return x;
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
-zero = 0.0,
-one = 1.0,
-two = 2.0,
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
huge = 1.0e300,
tiny = 1.0e-300,
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
huge = 1.0e300,
tiny = 1.0e-300,
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
- /* y == zero: x**0 = 1 */
+ /* y == 0.0: x**0 = 1 */
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3ff00000 && lx == 0)
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3ff00000 && lx == 0)
- /* y != zero: result is NaN if either arg is NaN */
+ /* y != 0.0: result is NaN if either arg is NaN */
if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) ||
iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0))
return (x+0.0) + (y+0.0);
if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) ||
iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0))
return (x+0.0) + (y+0.0);
if (ly == 0) {
if (iy == 0x7ff00000) { /* y is +-inf */
if (((ix-0x3ff00000)|lx) == 0) /* (-1)**+-inf is 1 */
if (ly == 0) {
if (iy == 0x7ff00000) { /* y is +-inf */
if (((ix-0x3ff00000)|lx) == 0) /* (-1)**+-inf is 1 */
else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
- return hy >= 0 ? y : zero;
+ return hy >= 0 ? y : 0.0;
else /* (|x|<1)**+-inf = 0,inf */
else /* (|x|<1)**+-inf = 0,inf */
- return hy < 0 ? -y : zero;
+ return hy < 0 ? -y : 0.0;
}
if (iy == 0x3ff00000) { /* y is +-1 */
if (hy < 0)
}
if (iy == 0x3ff00000) { /* y is +-1 */
if (hy < 0)
return x;
}
if (hy == 0x40000000) /* y is 2 */
return x;
}
if (hy == 0x40000000) /* y is 2 */
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */
z = ax;
if (hy < 0) /* z = (1/|x|) */
if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */
z = ax;
if (hy < 0) /* z = (1/|x|) */
if (hx < 0) {
if (((ix-0x3ff00000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
if (hx < 0) {
if (((ix-0x3ff00000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
if ((n|yisint) == 0)
return (x-x)/(x-x);
if ((n|yisint) == 0)
return (x-x)/(x-x);
- s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+ s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */
- s = -one;/* (-ve)**(odd int) */
+ s = -1.0;/* (-ve)**(odd int) */
/* |y| is huge */
if (iy > 0x41e00000) { /* if |y| > 2**31 */
/* |y| is huge */
if (iy > 0x41e00000) { /* if |y| > 2**31 */
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 */
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 */
+ t = ax - 1.0; /* 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;
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;
/* 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 */
/* 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 */
ss = u*v;
s_h = ss;
SET_LOW_WORD(s_h, 0);
/* t_h=ax+bp[k] High */
ss = u*v;
s_h = ss;
SET_LOW_WORD(s_h, 0);
/* t_h=ax+bp[k] High */
SET_HIGH_WORD(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);
SET_HIGH_WORD(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);
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 */
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 */
SET_HIGH_WORD(t, n & ~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if (j < 0)
SET_HIGH_WORD(t, n & ~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if (j < 0)
w = v - (z-u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
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);
+ r = (z*t1)/(t1-2.0) - (w + z*w);
+ z = 1.0 - (r-z);
GET_HIGH_WORD(j, z);
j += n<<20;
if ((j>>20) <= 0) /* subnormal output */
GET_HIGH_WORD(j, z);
j += n<<20;
if ((j>>20) <= 0) /* subnormal output */
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
bp[] = {1.0, 1.5,},
dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
-zero = 0.0,
-one = 1.0,
-two = 2.0,
two24 = 16777216.0, /* 0x4b800000 */
huge = 1.0e30,
tiny = 1.0e-30,
two24 = 16777216.0, /* 0x4b800000 */
huge = 1.0e30,
tiny = 1.0e-30,
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
- /* y == zero: x**0 = 1 */
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3f800000)
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3f800000)
- /* y != zero: result is NaN if either arg is NaN */
+ /* y != 0: result is NaN if either arg is NaN */
if (ix > 0x7f800000 || iy > 0x7f800000)
return (x+0.0f) + (y+0.0f);
if (ix > 0x7f800000 || iy > 0x7f800000)
return (x+0.0f) + (y+0.0f);
/* special value of y */
if (iy == 0x7f800000) { /* y is +-inf */
if (ix == 0x3f800000) /* (-1)**+-inf is 1 */
/* special value of y */
if (iy == 0x7f800000) { /* y is +-inf */
if (ix == 0x3f800000) /* (-1)**+-inf is 1 */
else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
- return hy >= 0 ? y : zero;
+ return hy >= 0 ? y : 0.0f;
else /* (|x|<1)**+-inf = 0,inf */
else /* (|x|<1)**+-inf = 0,inf */
- return hy < 0 ? -y : zero;
+ return hy < 0 ? -y : 0.0f;
}
if (iy == 0x3f800000) { /* y is +-1 */
if (hy < 0)
}
if (iy == 0x3f800000) { /* y is +-1 */
if (hy < 0)
return x;
}
if (hy == 0x40000000) /* y is 2 */
return x;
}
if (hy == 0x40000000) /* y is 2 */
if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
z = ax;
if (hy < 0) /* z = (1/|x|) */
if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
z = ax;
if (hy < 0) /* z = (1/|x|) */
if (hx < 0) {
if (((ix-0x3f800000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
if (hx < 0) {
if (((ix-0x3f800000)|yisint) == 0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
if ((n|yisint) == 0)
return (x-x)/(x-x);
if ((n|yisint) == 0)
return (x-x)/(x-x);
- sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+ sn = 1.0f; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */
if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */
/* |y| is huge */
if (iy > 0x4d000000) { /* if |y| > 2**27 */
/* |y| is huge */
if (iy > 0x4d000000) { /* if |y| > 2**27 */
/* compute s = 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 */
/* compute s = 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 */
s = u*v;
s_h = s;
GET_FLOAT_WORD(is, s_h);
s = u*v;
s_h = s;
GET_FLOAT_WORD(is, s_h);
w = v - (z - u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
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);
+ r = (z*t1)/(t1-2.0f) - (w+z*w);
+ z = 1.0f - (r - z);
GET_FLOAT_WORD(j, z);
j += n<<23;
if ((j>>23) <= 0) /* subnormal output */
GET_FLOAT_WORD(j, z);
j += n<<23;
if ((j>>23) <= 0) /* subnormal output */
#define MANL_SHIFT (LDBL_MANL_SIZE - 1)
#define MANL_SHIFT (LDBL_MANL_SIZE - 1)
-static const long double Zero[] = {0.0L, -0.0L};
+static const long double Zero[] = {0.0, -0.0};
/*
* Return the IEEE remainder and set *quo to the last n bits of the
/*
* Return the IEEE remainder and set *quo to the last n bits of the
* If the result is +-0, then it must have the same sign as x, but
* the above calculation doesn't always give this. Fix up the sign.
*/
* If the result is +-0, then it must have the same sign as x, but
* the above calculation doesn't always give this. Fix up the sign.
*/
- if (ex < BIAS && x == 0.0L)
+ if (ex < BIAS && x == 0.0)
return zero[sign];
return x;
return zero[sign];
return x;
return x * 0x1p-16382L;
}
}
return x * 0x1p-16382L;
}
}
scale.bits.exp = 0x3fff + n;
return x * scale.e;
}
scale.bits.exp = 0x3fff + n;
return x * scale.e;
}
-static const double one = 1.0, huge = 1.0e307;
+static const double huge = 1.0e307;
if (ix < 0x40360000) { /* |x|<22 */
if (ix < 0x3e300000) /* |x|<2**-28 */
/* raise inexact, return x */
if (ix < 0x40360000) { /* |x|<22 */
if (ix < 0x3e300000) /* |x|<2**-28 */
/* raise inexact, return x */
return x;
t = expm1(fabs(x));
if (ix < 0x3ff00000)
return x;
t = expm1(fabs(x));
if (ix < 0x3ff00000)
- return h*(2.0*t - t*t/(t+one));
- return h*(t + t/(t+one));
+ return h*(2.0*t - t*t/(t+1.0));
+ return h*(t + t/(t+1.0));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
-static const float one = 1.0, huge = 1.0e37;
+static const float huge = 1.0e37;
if (ix < 0x41100000) { /* |x|<9 */
if (ix < 0x39800000) /* |x|<2**-12 */
/* raise inexact, return x */
if (ix < 0x41100000) { /* |x|<9 */
if (ix < 0x39800000) /* |x|<2**-12 */
/* raise inexact, return x */
return x;
t = expm1f(fabsf(x));
if (ix < 0x3f800000)
return x;
t = expm1f(fabsf(x));
if (ix < 0x3f800000)
- return h*(2.0f*t - t*t/(t+one));
- return h*(t + t/(t+one));
+ return h*(2.0f*t - t*t/(t+1.0f));
+ return h*(t + t/(t+1.0f));
}
/* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */
}
/* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */
return sinh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
return sinh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double one = 1.0, huge = 1.0e4931L;
+static const long double huge = 1.0e4931L;
long double sinhl(long double x)
{
long double sinhl(long double x)
{
/* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */
if (ix < 0x3fdf) /* |x|<2**-32 */
/* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x| < 25 */
if (ix < 0x3fdf) /* |x|<2**-32 */
return x;/* sinh(tiny) = tiny with inexact */
t = expm1l(fabsl(x));
if (ix < 0x3fff)
return x;/* sinh(tiny) = tiny with inexact */
t = expm1l(fabsl(x));
if (ix < 0x3fff)
- return h*(2.0*t - t*t/(t + one));
- return h*(t + t/(t + one));
+ return h*(2.0*t - t*t/(t + 1.0));
+ return h*(t + t/(t + 1.0));
}
/* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
}
/* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
-static const double one = 1.0, tiny = 1.0e-300;
+static const double tiny = 1.0e-300;
/* use floating add to find out rounding direction */
if ((ix0|ix1) != 0) {
/* use floating add to find out rounding direction */
if ((ix0|ix1) != 0) {
- z = one - tiny; /* raise inexact flag */
- if (z >= one) {
- z = one + tiny;
+ z = 1.0 - tiny; /* raise inexact flag */
+ if (z >= 1.0) {
+ z = 1.0 + tiny;
if (q1 == (uint32_t)0xffffffff) {
q1 = 0;
q++;
if (q1 == (uint32_t)0xffffffff) {
q1 = 0;
q++;
if (q1 == (uint32_t)0xfffffffe)
q++;
q1 += 2;
if (q1 == (uint32_t)0xfffffffe)
q++;
q1 += 2;
-static const float one = 1.0, tiny = 1.0e-30;
+static const float tiny = 1.0e-30;
/* use floating add to find out rounding direction */
if (ix != 0) {
/* use floating add to find out rounding direction */
if (ix != 0) {
- z = one - tiny; /* raise inexact flag */
- if (z >= one) {
- z = one + tiny;
- if (z > one)
+ z = 1.0f - tiny; /* raise inexact flag */
+ if (z >= 1.0f) {
+ z = 1.0f + tiny;
+ if (z > 1.0f)
-static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
+static const double tiny = 1.0e-300, huge = 1.0e300;
/* x is INF or NaN */
if (ix >= 0x7ff00000) {
if (jx >= 0)
/* x is INF or NaN */
if (ix >= 0x7ff00000) {
if (jx >= 0)
- return one/x + one; /* tanh(+-inf)=+-1 */
+ return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
- return one/x - one; /* tanh(NaN) = NaN */
+ return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
}
if (ix < 0x40360000) { /* |x| < 22 */
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* tanh(tiny) = tiny with inexact */
}
if (ix < 0x40360000) { /* |x| < 22 */
if (ix < 0x3e300000) { /* |x| < 2**-28 */
/* tanh(tiny) = tiny with inexact */
return x;
}
if (ix >= 0x3ff00000) { /* |x| >= 1 */
return x;
}
if (ix >= 0x3ff00000) { /* |x| >= 1 */
- t = expm1(two*fabs(x));
- z = one - two/(t+two);
+ t = expm1(2.0f*fabs(x));
+ z = 1.0f - 2.0f/(t+2.0f);
- t = expm1(-two*fabs(x));
- z= -t/(t+two);
+ t = expm1(-2.0f*fabs(x));
+ z= -t/(t+2.0f);
}
} else { /* |x| >= 22, return +-1 */
}
} else { /* |x| >= 22, return +-1 */
- z = one - tiny; /* raise inexact */
+ z = 1.0f - tiny; /* raise inexact */
}
return jx >= 0 ? z : -z;
}
}
return jx >= 0 ? z : -z;
}
-static const float one = 1.0, two = 2.0, tiny = 1.0e-30, huge = 1.0e30;
+static const float
+tiny = 1.0e-30,
+huge = 1.0e30;
/* x is INF or NaN */
if(ix >= 0x7f800000) {
if (jx >= 0)
/* x is INF or NaN */
if(ix >= 0x7f800000) {
if (jx >= 0)
- return one/x + one; /* tanh(+-inf)=+-1 */
+ return 1.0f/x + 1.0f; /* tanh(+-inf)=+-1 */
- return one/x - one; /* tanh(NaN) = NaN */
+ return 1.0f/x - 1.0f; /* tanh(NaN) = NaN */
}
if (ix < 0x41100000) { /* |x| < 9 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* tanh(tiny) = tiny with inexact */
}
if (ix < 0x41100000) { /* |x| < 9 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
/* tanh(tiny) = tiny with inexact */
return x;
}
if (ix >= 0x3f800000) { /* |x|>=1 */
return x;
}
if (ix >= 0x3f800000) { /* |x|>=1 */
- t = expm1f(two*fabsf(x));
- z = one - two/(t+two);
+ t = expm1f(2.0f*fabsf(x));
+ z = 1.0f - 2.0f/(t+2.0f);
- t = expm1f(-two*fabsf(x));
- z = -t/(t+two);
+ t = expm1f(-2.0f*fabsf(x));
+ z = -t/(t+2.0f);
}
} else { /* |x| >= 9, return +-1 */
}
} else { /* |x| >= 9, return +-1 */
- z = one - tiny; /* raise inexact */
+ z = 1.0f - tiny; /* raise inexact */
}
return jx >= 0 ? z : -z;
}
}
return jx >= 0 ? z : -z;
}
return tanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
return tanh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
+static const long double tiny = 1.0e-4900L;
long double tanhl(long double x)
{
long double tanhl(long double x)
{
if (ix == 0x7fff) {
/* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (se & 0x8000)
if (ix == 0x7fff) {
/* for NaN it's not important which branch: tanhl(NaN) = NaN */
if (se & 0x8000)
- return one/x-one; /* tanhl(-inf)= -1; */
- return one/x+one; /* tanhl(+inf)= +1 */
+ return 1.0/x-1.0; /* tanhl(-inf)= -1; */
+ return 1.0/x+1.0; /* tanhl(+inf)= +1 */
if ((ix|jj0|jj1) == 0) /* x == +- 0 */
return x;
if (ix < 0x3fc8) /* |x| < 2**-55 */
if ((ix|jj0|jj1) == 0) /* x == +- 0 */
return x;
if (ix < 0x3fc8) /* |x| < 2**-55 */
- return x*(one+tiny); /* tanh(small) = small */
+ return x*(1.0+tiny); /* tanh(small) = small */
if (ix >= 0x3fff) { /* |x| >= 1 */
if (ix >= 0x3fff) { /* |x| >= 1 */
- t = expm1l(two*fabsl(x));
- z = one - two/(t+two);
+ t = expm1l(2.0*fabsl(x));
+ z = 1.0 - 2.0/(t+2.0);
- t = expm1l(-two*fabsl(x));
- z = -t/(t+two);
+ t = expm1l(-2.0*fabsl(x));
+ z = -t/(t+2.0);
}
/* |x| > 23, return +-1 */
} else {
}
/* |x| > 23, return +-1 */
} else {
- z = one - tiny; /* raise inexact flag */
+ z = 1.0 - tiny; /* raise inexact flag */
}
return se & 0x8000 ? -z : z;
}
}
return se & 0x8000 ? -z : z;
}
/* For large x, use rational coefficients from the analytical expansion. */
/* For large x, use rational coefficients from the analytical expansion. */
w = (((((6.97281375836585777429E-5L * w
+ 7.84039221720066627474E-4L) * w
- 2.29472093621399176955E-4L) * w
- 2.68132716049382716049E-3L) * w
+ 3.47222222222222222222E-3L) * w
+ 8.33333333333333333333E-2L) * w
w = (((((6.97281375836585777429E-5L * w
+ 7.84039221720066627474E-4L) * w
- 2.29472093621399176955E-4L) * w
- 2.68132716049382716049E-3L) * w
+ 3.47222222222222222222E-3L) * w
+ 8.33333333333333333333E-2L) * w
- w = 1.0L + w * __polevll(w, STIR, 8);
+ w = 1.0 + w * __polevll(w, STIR, 8);
y = expl(x);
if (x > MAXSTIR) { /* Avoid overflow in pow() */
v = powl(x, 0.5L * x - 0.25L);
y = expl(x);
if (x > MAXSTIR) { /* Avoid overflow in pow() */
v = powl(x, 0.5L * x - 0.25L);
if (x == -INFINITY)
return x - x;
q = fabsl(x);
if (x == -INFINITY)
return x - x;
q = fabsl(x);
if (q > MAXGAML)
goto goverf;
if (q > MAXGAML)
goto goverf;
p = floorl(q);
if (p == q)
return (x - x) / (x - x);
p = floorl(q);
if (p == q)
return (x - x) / (x - x);
signgam = -1;
z = q - p;
if (z > 0.5L) {
signgam = -1;
z = q - p;
if (z > 0.5L) {
z = q - p;
}
z = q * sinl(PIL * z);
z = q - p;
}
z = q * sinl(PIL * z);
- z = 1.0L;
- while (x >= 3.0L) {
- x -= 1.0L;
+ z = 1.0;
+ while (x >= 3.0) {
+ x -= 1.0;
z *= x;
}
while (x < -0.03125L) {
z /= x;
z *= x;
}
while (x < -0.03125L) {
z /= x;
}
if (x <= 0.03125L)
goto small;
}
if (x <= 0.03125L)
goto small;
p = __polevll(x, P, 7);
q = __polevll(x, Q, 8);
z = z * p / q;
p = __polevll(x, P, 7);
q = __polevll(x, Q, 8);
z = z * p / q;
return (x - x) / (x - x);
return (x - x) / (x - x);
x = -x;
q = z / (x * __polevll(x, SN, 8));
signgam = -1;
x = -x;
q = z / (x * __polevll(x, SN, 8));
signgam = -1;
projects / musl / commitdiff