Use OpenBSD 4.9's libm

This version still has all the double function (and not the long-double
ones). And apart from this, it's virtually identical to the 6.0 code.
This commit is contained in:
Romain Goyet
2017-11-08 16:14:24 +01:00
parent 247ea217fa
commit 3adc5d9fac
45 changed files with 141 additions and 367 deletions

View File

@@ -93,7 +93,6 @@ objs += $(addprefix liba/src/external/openbsd/, \
s_fabs.o \
s_floor.o \
s_log1p.o \
s_nan.o \
s_round.o \
s_scalbn.o \
s_sin.o \
@@ -103,36 +102,7 @@ objs += $(addprefix liba/src/external/openbsd/, \
)
liba/src/external/openbsd/%.o: SFLAGS := -Iliba/src/external/openbsd/include $(SFLAGS)
liba/src/external/openbsd/e_lgammaf_r.o: CFLAGS += -w
liba/src/external/openbsd/s_log1pf.o: CFLAGS += -w
liba/src/external/openbsd/s_scalbnf.o: CFLAGS += -w
liba/src/external/openbsd/e_acosh.o: CFLAGS += -w
liba/src/external/openbsd/e_atanh.o: CFLAGS += -w
liba/src/external/openbsd/e_cosh.o: CFLAGS += -w
liba/src/external/openbsd/e_exp.o: CFLAGS += -w
liba/src/external/openbsd/e_log.o: CFLAGS += -w
liba/src/external/openbsd/e_lgamma_r.o: CFLAGS += -w
liba/src/external/openbsd/e_pow.o: CFLAGS += -w
liba/src/external/openbsd/e_rem_pio2.o: CFLAGS += -w
liba/src/external/openbsd/e_sinh.o: CFLAGS += -w
liba/src/external/openbsd/k_rem_pio2.o: CFLAGS += -w
liba/src/external/openbsd/k_rem_pio2f.o: CFLAGS += -w
liba/src/external/openbsd/s_asinh.o: CFLAGS += -w
liba/src/external/openbsd/s_log1p.o: CFLAGS += -w
liba/src/external/openbsd/s_scalbn.o: CFLAGS += -w
liba/src/external/openbsd/s_tanh.o: CFLAGS += -w
liba/src/external/openbsd/w_lgamma.o: CFLAGS += -w
# some openbsd classes are throwing implicit declaration warnings
ifeq ($(DEBUG),1)
# OpenBSD uses double constants ("0.5" instead of "0.5f") in single-precision
# code. That's annoying because Clang rightfully decides to emit double-to-float
# aeabi conversions when building in -O0 mode, and we really don't want to code
# such functions. A simple workaround is to always build those files -Os.
liba/src/external/openbsd/e_expf.o: CFLAGS += -Os
liba/src/external/openbsd/s_expm1f.o: CFLAGS += -Os
liba/src/external/openbsd/s_log1pf.o: CFLAGS += -Os
liba/src/external/openbsd/s_roundf.o: CFLAGS += -Os
endif
liba/src/external/openbsd/%.o: CFLAGS += -w
tests += $(addprefix liba/test/, \
aeabi.c \

View File

@@ -1,7 +1,7 @@
"softfloat": Unaltered files from SoftFloat v3c.
http://www.jhauser.us/arithmetic/SoftFloat.html
"openbsd": Unaltered files from OpenBSD 6.0. Original path is "lib/libm/src"
"openbsd": Unaltered files from OpenBSD 4.9. Original path is "lib/libm/src"
"openbsd/include": Compatibility headers needed to build files from OpenBSD.
http://www.openbsd.org
@@ -11,5 +11,9 @@ single-precision variant (which fdlibm itself lacks). And more interestingly,
this variant does all its computation in single-precision (it never upgrades to
double-precision). In our case, this is very interesting because we're doing
single-precision computation when we better performance at the cost of accuracy,
and in our cse our hardware has a single-precision FPU so switching to double
and in our case our hardware has a single-precision FPU so switching to double
can yield much lower performances.
Last but not least, we're using OpenBSD 4.9 and not the latest one because this
version still has all the double precision code (some of it got removed later on
in favor of a "long double" version).

View File

@@ -34,6 +34,7 @@
* Function needed: sqrt
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -100,6 +101,8 @@ acos(double x)
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(acosl, acos);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(acosl, acos);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -24,9 +24,7 @@
* acosh(NaN) is NaN without signal.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -57,7 +55,3 @@ acosh(double x)
return log1p(t+sqrt(2.0*t+t*t));
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(acoshl, acosh);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -40,6 +40,7 @@
*
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -80,12 +81,12 @@ asin(double x)
} else if (ix<0x3fe00000) { /* |x|<0.5 */
if(ix<0x3e400000) { /* if |x| < 2**-27 */
if(huge+x>one) return x;/* return x with inexact if x!=0*/
}
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x+x*w;
} else
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x+x*w;
}
/* 1> |x|>= 0.5 */
w = one-fabs(x);
@@ -108,6 +109,8 @@ asin(double x)
if(hx>0) return t; else return -t;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(asinl, asin);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(asinl, asin);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -49,12 +49,12 @@ asinf(float x)
} else if (ix<0x3f000000) { /* |x|<0.5 */
if(ix<0x32000000) { /* if |x| < 2**-27 */
if(huge+x>one) 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)));
w = p/q;
return x+x*w;
} else
t = x*x;
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
w = p/q;
return x+x*w;
}
/* 1> |x|>= 0.5 */
w = one-fabsf(x);

View File

@@ -28,9 +28,7 @@
*
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double one = 1.0, huge = 1e300;
@@ -57,7 +55,3 @@ atanh(double x)
t = 0.5*log1p((x+x)/(one-x));
if(hx>=0) return t; else return -t;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(atanhl, atanh);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -31,9 +31,7 @@
* only cosh(0)=1 is exact for finite x.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double one = 1.0, half=0.5, huge = 1.0e300;
@@ -81,7 +79,3 @@ cosh(double x)
/* |x| > overflowthresold, cosh(x) overflow */
return huge*huge;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(coshl, cosh);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -73,9 +73,7 @@
* to produce the hexadecimal values shown.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -155,7 +153,3 @@ exp(double x) /* default IEEE double exp */
return y*twom1000;
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(expl, exp);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -45,7 +45,6 @@ expf(float x) /* default IEEE double exp */
GET_FLOAT_WORD(hx,x);
xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
k = 0;
/* filter out non-finite argument */
if(hx >= 0x42b17218) { /* if |x|>=88.721... */
@@ -72,6 +71,7 @@ expf(float x) /* default IEEE double exp */
else if(hx < 0x31800000) { /* when |x|<2**-28 */
if(huge+x>one) return one+x;/* trigger inexact */
}
else k = 0;
/* x is now in primary range */
t = x*x;

View File

@@ -208,11 +208,7 @@ lgamma_r(double x, int *signgamp)
*signgamp = 1;
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x*x;
if((ix|lx)==0) {
if(hx<0)
*signgamp = -1;
return one/zero;
}
if((ix|lx)==0) return one/zero;
if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */
if(hx<0) {
*signgamp = -1;

View File

@@ -144,11 +144,7 @@ lgammaf_r(float x, int *signgamp)
*signgamp = 1;
ix = hx&0x7fffffff;
if(ix>=0x7f800000) return x*x;
if(ix==0) {
if(hx<0)
*signgamp = -1;
return one/zero;
}
if(ix==0) return one/zero;
if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */
if(hx<0) {
*signgamp = -1;

View File

@@ -11,7 +11,7 @@
*/
/* log(x)
* Return the logarithm of x
* Return the logrithm of x
*
* Method :
* 1. Argument Reduction: find k and f such that
@@ -61,9 +61,7 @@
* to produce the hexadecimal values shown.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -130,7 +128,3 @@ log(double x)
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(logl, log);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -43,9 +43,7 @@
* shown.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -82,7 +80,3 @@ log10(double x)
z = y*log10_2lo + ivln10*log(x);
return z+y*log10_2hi;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(log10l, log10);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -54,8 +54,8 @@ logf(float x)
k += (i>>23);
f = x-(float)1.0;
if((0x007fffff&(15+ix))<16) { /* |f| < 2**-20 */
if(f==zero) {if(k==0) return zero; else {dk=(float)k;
return dk*ln2_hi+dk*ln2_lo;}}
if(f==zero) if(k==0) return zero; else {dk=(float)k;
return dk*ln2_hi+dk*ln2_lo;}
R = f*f*((float)0.5-(float)0.33333333333333333*f);
if(k==0) return f-R; else {dk=(float)k;
return dk*ln2_hi-((R-dk*ln2_lo)-f);}

View File

@@ -24,13 +24,13 @@
* Special cases:
* 1. (anything) ** 0 is 1
* 2. (anything) ** 1 is itself
* 3. (anything except 1) ** NAN is NAN
* 3. (anything) ** NAN is NAN
* 4. NAN ** (anything except 0) is NAN
* 5. +-(|x| > 1) ** +INF is +INF
* 6. +-(|x| > 1) ** -INF is +0
* 7. +-(|x| < 1) ** +INF is +0
* 8. +-(|x| < 1) ** -INF is +INF
* 9. +-1 ** +-INF is 1
* 9. +-1 ** +-INF is NAN
* 10. +0 ** (+anything except 0, NAN) is +0
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF
@@ -55,9 +55,7 @@
* to produce the hexadecimal values shown.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -109,9 +107,6 @@ pow(double x, double y)
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
/* x==1: 1**y = 1, even if y is NaN */
if (hx==0x3ff00000 && lx == 0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
@@ -141,7 +136,7 @@ pow(double x, double y)
if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if(((ix-0x3ff00000)|lx)==0)
return one; /* (-1)**+-inf is 1 */
return y - y; /* inf**+-1 is NaN */
else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */
@@ -300,7 +295,3 @@ pow(double x, double y)
else SET_HIGH_WORD(z,j);
return s*z;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(powl, pow);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -64,9 +64,6 @@ powf(float x, float y)
/* y==zero: x**0 = 1 */
if(iy==0) return one;
/* x==1: 1**y = 1, even if y is NaN */
if (hx==0x3f800000) return one;
/* +-NaN return x+y */
if(ix > 0x7f800000 ||
iy > 0x7f800000)
@@ -90,7 +87,7 @@ powf(float x, float y)
/* special value of y */
if (iy==0x7f800000) { /* y is +-inf */
if (ix==0x3f800000)
return one; /* (-1)**+-inf is NaN */
return y - y; /* inf**+-1 is NaN */
else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */

View File

@@ -28,9 +28,7 @@
* only sinh(0)=0 is exact for finite x.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double one = 1.0, shuge = 1.0e307;
@@ -74,7 +72,3 @@ sinh(double x)
/* |x| > overflowthresold, sinh(x) overflow */
return x*shuge;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(sinhl, sinh);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -80,6 +80,7 @@
*---------------
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -441,6 +442,8 @@ B. sqrt(x) by Reciproot Iteration
*/
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(sqrtl, sqrt);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(sqrtl, sqrt);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -1,5 +1,10 @@
#include_next <math.h>
/* s_roundf.c uses isinff and isnanf. Let's correct this */
#define isinff(x) __builtin_isinf(x)
#define isnanf(x) __builtin_isnan(x)
/* In accordance with the C99 standard, we've defined libm function as macros to
* leverage compiler optimizations. When comes the time to actually implement
* those functions, we don't want the macro to be active. OpenBSD doesn't use

View File

@@ -0,0 +1 @@
// OpenBSD requires this file to be present

View File

@@ -21,9 +21,7 @@
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -53,7 +51,3 @@ asinh(double x)
}
if(hx>0) return w; else return -w;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(asinhl, asinh);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -30,6 +30,7 @@
* to produce the hexadecimal values shown.
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -116,6 +117,8 @@ atan(double x)
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(atanl, atan);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(atanl, atan);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -19,9 +19,7 @@
* Inexact flag raised if x not equal to ceil(x).
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double huge = 1.0e300;
@@ -68,7 +66,3 @@ ceil(double x)
INSERT_WORDS(x,i0,i1);
return x;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(ceill, ceil);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -46,4 +46,4 @@ ceilf(float x)
}
SET_FLOAT_WORD(x,i0);
return x;
}
}

View File

@@ -16,6 +16,7 @@
* with the sign bit of y.
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -31,6 +32,8 @@ copysign(double x, double y)
return x;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(copysignl, copysign);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(copysignl, copysign);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -41,6 +41,7 @@
* TRIG(x) returns trig(x) nearly rounded
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -75,6 +76,8 @@ cos(double x)
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(cosl, cos);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(cosl, cos);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -105,9 +105,7 @@
* to produce the hexadecimal values shown.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -188,10 +186,9 @@ expm1(double x)
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return 0.5*(x-e)-0.5;
if(k==1) {
if(k==1)
if(x < -0.25) return -2.0*(e-(x+0.5));
else return one+2.0*(x-e);
}
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
u_int32_t high;
y = one-(e-x);
@@ -217,7 +214,3 @@ expm1(double x)
}
return y;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(expm1l, expm1);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -91,10 +91,9 @@ expm1f(float x)
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return (float)0.5*(x-e)-(float)0.5;
if(k==1) {
if(k==1)
if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5));
else return one+(float)2.0*(x-e);
}
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
int32_t i;
y = one-(e-x);

View File

@@ -5,7 +5,7 @@
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -14,7 +14,10 @@
* fabs(x) returns the absolute value of x.
*/
#include "math.h"
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
double
@@ -25,3 +28,9 @@ fabs(double x)
SET_HIGH_WORD(x,high&0x7fffffff);
return x;
}
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(fabsl, fabs);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -19,9 +19,7 @@
* Inexact flag raised if x not equal to floor(x).
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double huge = 1.0e300;
@@ -69,7 +67,3 @@ floor(double x)
INSERT_WORDS(x,i0,i1);
return x;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(floorl, floor);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -35,6 +35,6 @@ frexpf(float x, int *eptr)
}
*eptr += (ix>>23)-126;
hx = (hx&0x807fffff)|0x3f000000;
SET_FLOAT_WORD(x,hx);
*(int*)&x = hx;
return x;
}

View File

@@ -75,9 +75,7 @@
* See HP-15C Advanced Functions Handbook, p.193.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double
@@ -157,7 +155,3 @@ log1p(double x)
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(log1pl, log1p);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -1,127 +0,0 @@
/* $OpenBSD: s_nan.c,v 1.13 2014/07/21 01:51:11 guenther Exp $ */
/*-
* Copyright (c) 2007 David Schultz
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* $FreeBSD: src/lib/msun/src/s_nan.c,v 1.2 2007/12/18 23:46:32 das Exp $
*/
#include <sys/types.h>
#include <ctype.h>
#include <endian.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <strings.h>
#include "math_private.h"
/*
* OpenBSD's ctype doesn't have digittoint, so we define it here.
*/
static int
_digittoint(int c)
{
if (!isxdigit(c))
return 0;
if (isdigit(c))
return c - '0';
else
return isupper(c) ? c - 'A' + 10 : c - 'a' + 10;
}
/*
* Scan a string of hexadecimal digits (the format nan(3) expects) and
* make a bit array (using the local endianness). We stop when we
* encounter an invalid character, NUL, etc. If we overflow, we do
* the same as gcc's __builtin_nan(), namely, discard the high order bits.
*
* The format this routine accepts needs to be compatible with what is used
* in contrib/gdtoa/hexnan.c (for strtod/scanf) and what is used in
* __builtin_nan(). In fact, we're only 100% compatible for strings we
* consider valid, so we might be violating the C standard. But it's
* impossible to use nan(3) portably anyway, so this seems good enough.
*/
void
_scan_nan(uint32_t *words, int num_words, const char *s)
{
int si; /* index into s */
int bitpos; /* index into words (in bits) */
bzero(words, num_words * sizeof(uint32_t));
/* Allow a leading '0x'. (It's expected, but redundant.) */
if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
s += 2;
/* Scan forwards in the string, looking for the end of the sequence. */
for (si = 0; isxdigit((unsigned char)s[si]); si++)
;
/* Scan backwards, filling in the bits in words[] as we go. */
#if BYTE_ORDER == LITTLE_ENDIAN
for (bitpos = 0; bitpos < 32 * num_words; bitpos += 4) {
#else
for (bitpos = 32 * num_words - 4; bitpos >= 0; bitpos -= 4) {
#endif
if (--si < 0)
break;
words[bitpos / 32] |= _digittoint((unsigned char)s[si]) << (bitpos % 32);
}
}
double
nan(const char *s)
{
union {
double d;
uint32_t bits[2];
} u;
_scan_nan(u.bits, 2, s);
#if BYTE_ORDER == LITTLE_ENDIAN
u.bits[1] |= 0x7ff80000;
#else
u.bits[0] |= 0x7ff80000;
#endif
return (u.d);
}
float
nanf(const char *s)
{
union {
float f;
uint32_t bits[1];
} u;
_scan_nan(u.bits, 1, s);
u.bits[0] |= 0x7fc00000;
return (u.f);
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(nanl, nan);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -1,4 +1,4 @@
/* $OpenBSD: s_round.c,v 1.6 2013/07/03 04:46:36 espie Exp $ */
/* $OpenBSD: s_round.c,v 1.1 2006/07/12 07:26:08 brad Exp $ */
/*-
* Copyright (c) 2003, Steven G. Kargl
@@ -26,9 +26,7 @@
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
double
@@ -51,7 +49,3 @@ round(double x)
return (-t);
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(roundl, round);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -34,7 +34,7 @@ roundf(float x)
{
float t;
if (isinf(x) || isnan(x))
if (isinff(x) || isnanf(x))
return (x);
if (x >= 0.0) {

View File

@@ -5,18 +5,19 @@
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
/*
* scalbn (double x, int n)
* scalbn(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* scalbn(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* exponentiation or a multiplication.
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -36,13 +37,13 @@ scalbn (double x, int n)
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
x *= two54;
GET_HIGH_WORD(hx,x);
k = ((hx&0x7ff00000)>>20) - 54;
k = ((hx&0x7ff00000)>>20) - 54;
if (n< -50000) return tiny*x; /*underflow*/
}
if (k==0x7ff) return x+x; /* NaN or Inf */
k = k+n;
k = k+n;
if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
if (k > 0) /* normal result */
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
@@ -55,9 +56,8 @@ scalbn (double x, int n)
return x*twom54;
}
double
ldexp(double x, int n)
{
return scalbn(x, n);
}
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(scalbnl, scalbn);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -1,4 +1,3 @@
/* $OpenBSD: s_signgam.c,v 1.2 2015/01/20 04:41:01 krw Exp $ */
#include "math.h"
#include "math_private.h"
int signgam = 0;

View File

@@ -41,6 +41,7 @@
* TRIG(x) returns trig(x) nearly rounded
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -75,6 +76,8 @@ sin(double x)
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(sinl, sin);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(sinl, sin);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -40,6 +40,7 @@
* TRIG(x) returns trig(x) nearly rounded
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
@@ -69,6 +70,8 @@ tan(double x)
}
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(tanl, tan);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(tanl, tan);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

View File

@@ -20,23 +20,21 @@
* x -x
* e + e
* 1. reduce x to non-negative by tanh(-x) = -tanh(x).
* 2. 0 <= x < 2**-55 : tanh(x) := x*(one+x)
* 2. 0 <= x <= 2**-55 : tanh(x) := x*(one+x)
* -t
* 2**-55 <= x < 1 : tanh(x) := -----; t = expm1(-2x)
* 2**-55 < x <= 1 : tanh(x) := -----; t = expm1(-2x)
* t + 2
* 2
* 1 <= x < 22.0 : tanh(x) := 1- ----- ; t=expm1(2x)
* 1 <= x <= 22.0 : tanh(x) := 1- ----- ; t=expm1(2x)
* t + 2
* 22.0 <= x <= INF : tanh(x) := 1.
* 22.0 < x <= INF : tanh(x) := 1.
*
* Special cases:
* tanh(NaN) is NaN;
* only tanh(0)=0 is exact for finite argument.
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
static const double one=1.0, two=2.0, tiny = 1.0e-300;
@@ -45,10 +43,10 @@ double
tanh(double x)
{
double t,z;
int32_t jx,ix,lx;
int32_t jx,ix;
/* High word of |x|. */
EXTRACT_WORDS(jx,lx,x);
GET_HIGH_WORD(jx,x);
ix = jx&0x7fffffff;
/* x is INF or NaN */
@@ -59,8 +57,6 @@ tanh(double x)
/* |x| < 22 */
if (ix < 0x40360000) { /* |x|<22 */
if ((ix | lx) == 0)
return x; /* x == +-0 */
if (ix<0x3c800000) /* |x|<2**-55 */
return x*(one+x); /* tanh(small) = small */
if (ix>=0x3ff00000) { /* |x|>=1 */
@@ -70,13 +66,9 @@ tanh(double x)
t = expm1(-two*fabs(x));
z= -t/(t+two);
}
/* |x| >= 22, return +-1 */
/* |x| > 22, return +-1 */
} else {
z = one - tiny; /* raised inexact flag */
}
return (jx>=0)? z: -z;
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(tanhl, tanh);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */

View File

@@ -35,8 +35,6 @@ tanhf(float x)
/* |x| < 22 */
if (ix < 0x41b00000) { /* |x|<22 */
if (ix == 0)
return x; /* x == +-0 */
if (ix<0x24000000) /* |x|<2**-55 */
return x*(one+x); /* tanh(small) = small */
if (ix>=0x3f800000) { /* |x|>=1 */

View File

@@ -10,6 +10,11 @@
* ====================================================
*/
#if 0
#include <sys/cdefs.h>
__FBSDID("$FreeBSD: src/lib/msun/src/s_truncf.c,v 1.1 2004/06/20 09:25:43 das Exp $");
#endif
/*
* truncf(x)
* Return x rounded toward 0 to integral value

View File

@@ -16,9 +16,7 @@
* Method: call lgamma_r
*/
#include <float.h>
#include <math.h>
#include "math.h"
#include "math_private.h"
extern int signgam;
@@ -28,7 +26,3 @@ lgamma(double x)
{
return lgamma_r(x,&signgam);
}
#if LDBL_MANT_DIG == DBL_MANT_DIG
__strong_alias(lgammal, lgamma);
#endif /* LDBL_MANT_DIG == DBL_MANT_DIG */