mirror of https://git.musl-libc.org/git/musl
Browse Source
modifications:
* avoid unsigned->signed integer conversion
* do not handle special cases when they work correctly anyway
* more strict threshold values (0x1p26 instead of 0x1p28 etc)
* smaller code, cleaner branching logic
* same precision as the old code:
acosh(x) has up to 2ulp error in [1,1.125]
asinh(x) has up to 1.6ulp error in [0.125,0.5], [-0.5,-0.125]
atanh(x) has up to 1.7ulp error in [0.125,0.5], [-0.5,-0.125]
rs-1.0
Szabolcs Nagy
14 years ago
9 changed files with 149 additions and 406 deletions
@ -1,54 +1,19 @@ |
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_acosh.c */ |
|||
/*
|
|||
* ==================================================== |
|||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|||
* |
|||
* Developed at SunSoft, a Sun Microsystems, Inc. business. |
|||
* Permission to use, copy, modify, and distribute this |
|||
* software is freely granted, provided that this notice |
|||
* is preserved. |
|||
* ==================================================== |
|||
* |
|||
*/ |
|||
/* acosh(x)
|
|||
* Method : |
|||
* Based on |
|||
* acosh(x) = log [ x + sqrt(x*x-1) ] |
|||
* we have |
|||
* acosh(x) := log(x)+ln2, if x is large; else |
|||
* acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else |
|||
* acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. |
|||
* |
|||
* Special cases: |
|||
* acosh(x) is NaN with signal if x<1. |
|||
* acosh(NaN) is NaN without signal. |
|||
*/ |
|||
|
|||
#include "libm.h" |
|||
|
|||
static const double |
|||
ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ |
|||
|
|||
/* acosh(x) = log(x + sqrt(x*x-1)) */ |
|||
double acosh(double x) |
|||
{ |
|||
double t; |
|||
int32_t hx; |
|||
uint32_t lx; |
|||
union {double f; uint64_t i;} u = {.f = x}; |
|||
unsigned e = u.i >> 52 & 0x7ff; |
|||
|
|||
/* x < 1 domain error is handled in the called functions */ |
|||
|
|||
EXTRACT_WORDS(hx, lx, x); |
|||
if (hx < 0x3ff00000) { /* x < 1 */ |
|||
return (x-x)/(x-x); |
|||
} else if (hx >= 0x41b00000) { /* x > 2**28 */ |
|||
if (hx >= 0x7ff00000) /* x is inf of NaN */ |
|||
return x+x; |
|||
return log(x) + ln2; /* acosh(huge) = log(2x) */ |
|||
} else if ((hx-0x3ff00000 | lx) == 0) { |
|||
return 0.0; /* acosh(1) = 0 */ |
|||
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */ |
|||
t = x*x; |
|||
return log(2.0*x - 1.0/(x+sqrt(t-1.0))); |
|||
} else { /* 1 < x < 2 */ |
|||
t = x-1.0; |
|||
return log1p(t + sqrt(2.0*t+t*t)); |
|||
} |
|||
if (e < 0x3ff + 1) |
|||
/* |x| < 2, up to 2ulp error in [1,1.125] */ |
|||
return log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1))); |
|||
if (e < 0x3ff + 26) |
|||
/* |x| < 0x1p26 */ |
|||
return log(2*x - 1/(x+sqrt(x*x-1))); |
|||
/* |x| >= 0x1p26 or nan */ |
|||
return log(x) + 0.693147180559945309417232121458176568; |
|||
} |
|||
|
|||
@ -1,42 +1,17 @@ |
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_acoshf.c */ |
|||
/*
|
|||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
|||
*/ |
|||
/*
|
|||
* ==================================================== |
|||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|||
* |
|||
* Developed at SunPro, a Sun Microsystems, Inc. business. |
|||
* Permission to use, copy, modify, and distribute this |
|||
* software is freely granted, provided that this notice |
|||
* is preserved. |
|||
* ==================================================== |
|||
*/ |
|||
|
|||
#include "libm.h" |
|||
|
|||
static const float |
|||
ln2 = 6.9314718246e-01; /* 0x3f317218 */ |
|||
|
|||
/* acosh(x) = log(x + sqrt(x*x-1)) */ |
|||
float acoshf(float x) |
|||
{ |
|||
float t; |
|||
int32_t hx; |
|||
union {float f; int32_t i;} u = {.f = x}; |
|||
|
|||
GET_FLOAT_WORD(hx, x); |
|||
if (hx < 0x3f800000) { /* x < 1 */ |
|||
return (x-x)/(x-x); |
|||
} else if (hx >= 0x4d800000) { /* x > 2**28 */ |
|||
if (hx >= 0x7f800000) /* x is inf of NaN */ |
|||
return x + x; |
|||
return logf(x) + ln2; /* acosh(huge)=log(2x) */ |
|||
} else if (hx == 0x3f800000) { |
|||
return 0.0f; /* acosh(1) = 0 */ |
|||
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */ |
|||
t = x*x; |
|||
return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f))); |
|||
} else { /* 1 < x < 2 */ |
|||
t = x-1.0f; |
|||
return log1pf(t + sqrtf(2.0f*t+t*t)); |
|||
} |
|||
if (u.i < 0x3f800000+(1<<23)) |
|||
/* x < 2, invalid if x < 1 or nan */ |
|||
/* up to 2ulp error in [1,1.125] */ |
|||
return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1))); |
|||
if (u.i < 0x3f800000+(12<<23)) |
|||
/* x < 0x1p12 */ |
|||
return logf(2*x - 1/(x+sqrtf(x*x-1))); |
|||
/* x >= 0x1p12 */ |
|||
return logf(x) + 0.693147180559945309417232121458176568f; |
|||
} |
|||
|
|||
@ -1,55 +1,28 @@ |
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_asinh.c */ |
|||
/*
|
|||
* ==================================================== |
|||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|||
* |
|||
* Developed at SunPro, a Sun Microsystems, Inc. business. |
|||
* Permission to use, copy, modify, and distribute this |
|||
* software is freely granted, provided that this notice |
|||
* is preserved. |
|||
* ==================================================== |
|||
*/ |
|||
/* asinh(x)
|
|||
* Method : |
|||
* Based on |
|||
* asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] |
|||
* we have |
|||
* asinh(x) := x if 1+x*x=1, |
|||
* := sign(x)*(log(x)+ln2)) for large |x|, else |
|||
* := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else |
|||
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) |
|||
*/ |
|||
|
|||
#include "libm.h" |
|||
|
|||
static const double |
|||
ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ |
|||
huge= 1.00000000000000000000e+300; |
|||
|
|||
/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ |
|||
double asinh(double x) |
|||
{ |
|||
double t,w; |
|||
int32_t hx,ix; |
|||
union {double f; uint64_t i;} u = {.f = x}; |
|||
unsigned e = u.i >> 52 & 0x7ff; |
|||
unsigned s = u.i >> 63; |
|||
|
|||
GET_HIGH_WORD(hx, x); |
|||
ix = hx & 0x7fffffff; |
|||
if (ix >= 0x7ff00000) /* x is inf or NaN */ |
|||
return x+x; |
|||
if (ix < 0x3e300000) { /* |x| < 2**-28 */ |
|||
/* return x inexact except 0 */ |
|||
if (huge+x > 1.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); |
|||
w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t)); |
|||
} else { /* 2.0 > |x| > 2**-28 */ |
|||
t = x*x; |
|||
w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t))); |
|||
/* |x| */ |
|||
u.i &= (uint64_t)-1/2; |
|||
x = u.f; |
|||
|
|||
if (e >= 0x3ff + 26) { |
|||
/* |x| >= 0x1p26 or inf or nan */ |
|||
x = log(x) + 0.693147180559945309417232121458176568; |
|||
} else if (e >= 0x3ff + 1) { |
|||
/* |x| >= 2 */ |
|||
x = log(2*x + 1/(sqrt(x*x+1)+x)); |
|||
} else if (e >= 0x3ff - 26) { |
|||
/* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */ |
|||
x = log1p(x + x*x/(sqrt(x*x+1)+1)); |
|||
} else { |
|||
/* |x| < 0x1p-26, raise inexact if x != 0 */ |
|||
FORCE_EVAL(x + 0x1p1000); |
|||
} |
|||
if (hx > 0) |
|||
return w; |
|||
return -w; |
|||
return s ? -x : x; |
|||
} |
|||
|
|||
@ -1,48 +1,28 @@ |
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_asinhf.c */ |
|||
/*
|
|||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
|||
*/ |
|||
/*
|
|||
* ==================================================== |
|||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|||
* |
|||
* Developed at SunPro, a Sun Microsystems, Inc. business. |
|||
* Permission to use, copy, modify, and distribute this |
|||
* software is freely granted, provided that this notice |
|||
* is preserved. |
|||
* ==================================================== |
|||
*/ |
|||
|
|||
#include "libm.h" |
|||
|
|||
static const float |
|||
ln2 = 6.9314718246e-01, /* 0x3f317218 */ |
|||
huge= 1.0000000000e+30; |
|||
|
|||
/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ |
|||
float asinhf(float x) |
|||
{ |
|||
float t,w; |
|||
int32_t hx,ix; |
|||
union {float f; uint32_t i;} u = {.f = x}; |
|||
uint32_t i = u.i & 0x7fffffff; |
|||
unsigned s = u.i >> 31; |
|||
|
|||
GET_FLOAT_WORD(hx, x); |
|||
ix = hx & 0x7fffffff; |
|||
if (ix >= 0x7f800000) /* x is inf or NaN */ |
|||
return x+x; |
|||
if (ix < 0x31800000) { /* |x| < 2**-28 */ |
|||
/* return x inexact except 0 */ |
|||
if (huge+x > 1.0f) |
|||
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 + 1.0f/(sqrtf(x*x+1.0f)+t)); |
|||
} else { /* 2.0 > |x| > 2**-28 */ |
|||
t = x*x; |
|||
w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t))); |
|||
/* |x| */ |
|||
u.i = i; |
|||
x = u.f; |
|||
|
|||
if (i >= 0x3f800000 + (12<<23)) { |
|||
/* |x| >= 0x1p12 or inf or nan */ |
|||
x = logf(x) + 0.693147180559945309417232121458176568f; |
|||
} else if (i >= 0x3f800000 + (1<<23)) { |
|||
/* |x| >= 2 */ |
|||
x = logf(2*x + 1/(sqrtf(x*x+1)+x)); |
|||
} else if (i >= 0x3f800000 - (12<<23)) { |
|||
/* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */ |
|||
x = log1pf(x + x*x/(sqrtf(x*x+1)+1)); |
|||
} else { |
|||
/* |x| < 0x1p-12, raise inexact if x!=0 */ |
|||
FORCE_EVAL(x + 0x1p120f); |
|||
} |
|||
if (hx > 0) |
|||
return w; |
|||
return -w; |
|||
return s ? -x : x; |
|||
} |
|||
|
|||
@ -1,58 +1,21 @@ |
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_atanh.c */ |
|||
/*
|
|||
* ==================================================== |
|||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|||
* |
|||
* Developed at SunSoft, a Sun Microsystems, Inc. business. |
|||
* Permission to use, copy, modify, and distribute this |
|||
* software is freely granted, provided that this notice |
|||
* is preserved. |
|||
* ==================================================== |
|||
* |
|||
*/ |
|||
/* atanh(x)
|
|||
* Method : |
|||
* 1.Reduced x to positive by atanh(-x) = -atanh(x) |
|||
* 2.For x>=0.5 |
|||
* 1 2x x |
|||
* atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) |
|||
* 2 1 - x 1 - x |
|||
* |
|||
* For x<0.5 |
|||
* atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) |
|||
* |
|||
* Special cases: |
|||
* atanh(x) is NaN if |x| > 1 with signal; |
|||
* atanh(NaN) is that NaN with no signal; |
|||
* atanh(+-1) is +-INF with signal. |
|||
* |
|||
*/ |
|||
|
|||
#include "libm.h" |
|||
|
|||
static const double huge = 1e300; |
|||
|
|||
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ |
|||
double atanh(double x) |
|||
{ |
|||
double t; |
|||
int32_t hx,ix; |
|||
uint32_t lx; |
|||
union {double f; uint64_t i;} u = {.f = x}; |
|||
unsigned e = u.i >> 52 & 0x7ff; |
|||
unsigned s = u.i >> 63; |
|||
|
|||
/* |x| */ |
|||
u.i &= (uint64_t)-1/2; |
|||
x = u.f; |
|||
|
|||
EXTRACT_WORDS(hx, lx, x); |
|||
ix = hx & 0x7fffffff; |
|||
if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */ |
|||
return (x-x)/(x-x); |
|||
if (ix == 0x3ff00000) |
|||
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; |
|||
t = 0.5*log1p(t + t*x/(1.0-x)); |
|||
} else |
|||
t = 0.5*log1p((x+x)/(1.0-x)); |
|||
if (hx >= 0) |
|||
return t; |
|||
return -t; |
|||
if (e < 0x3ff - 1) { |
|||
/* |x| < 0.5, up to 1.7ulp error */ |
|||
x = 0.5*log1p(2*x + 2*x*x/(1-x)); |
|||
} else { |
|||
x = 0.5*log1p(2*x/(1-x)); |
|||
} |
|||
return s ? -x : x; |
|||
} |
|||
|
|||
@ -1,42 +1,20 @@ |
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_atanhf.c */ |
|||
/*
|
|||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
|||
*/ |
|||
/*
|
|||
* ==================================================== |
|||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|||
* |
|||
* Developed at SunPro, a Sun Microsystems, Inc. business. |
|||
* Permission to use, copy, modify, and distribute this |
|||
* software is freely granted, provided that this notice |
|||
* is preserved. |
|||
* ==================================================== |
|||
*/ |
|||
|
|||
#include "libm.h" |
|||
|
|||
static const float huge = 1e30; |
|||
|
|||
/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ |
|||
float atanhf(float x) |
|||
{ |
|||
float t; |
|||
int32_t hx,ix; |
|||
union {float f; uint32_t i;} u = {.f = x}; |
|||
unsigned s = u.i >> 31; |
|||
|
|||
/* |x| */ |
|||
u.i &= 0x7fffffff; |
|||
x = u.f; |
|||
|
|||
GET_FLOAT_WORD(hx, x); |
|||
ix = hx & 0x7fffffff; |
|||
if (ix > 0x3f800000) /* |x| > 1 */ |
|||
return (x-x)/(x-x); |
|||
if (ix == 0x3f800000) |
|||
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; |
|||
t = 0.5f*log1pf(t + t*x/(1.0f-x)); |
|||
} else |
|||
t = 0.5f*log1pf((x+x)/(1.0f-x)); |
|||
if (hx >= 0) |
|||
return t; |
|||
return -t; |
|||
if (u.i < 0x3f800000 - (1<<23)) { |
|||
/* |x| < 0.5, up to 1.7ulp error */ |
|||
x = 0.5f*log1pf(2*x + 2*x*x/(1-x)); |
|||
} else { |
|||
x = 0.5f*log1pf(2*x/(1-x)); |
|||
} |
|||
return s ? -x : x; |
|||
} |
|||
|
|||
Loading…
Reference in new issue