about summary refs log tree commit diff
path: root/src/math
diff options
context:
space:
mode:
Diffstat (limited to 'src/math')
-rw-r--r--src/math/__rem_pio2.c14
-rw-r--r--src/math/__rem_pio2f.c13
-rw-r--r--src/math/__rem_pio2l.c6
-rw-r--r--src/math/ceil.c11
-rw-r--r--src/math/ceill.c12
-rw-r--r--src/math/floor.c11
-rw-r--r--src/math/floorl.c12
-rw-r--r--src/math/modfl.c10
-rw-r--r--src/math/rintl.c12
-rw-r--r--src/math/round.c11
-rw-r--r--src/math/roundf.c13
-rw-r--r--src/math/roundl.c12
-rw-r--r--src/math/truncl.c10
13 files changed, 89 insertions, 58 deletions
diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c
index 5fafc4a8..a40db9fc 100644
--- a/src/math/__rem_pio2.c
+++ b/src/math/__rem_pio2.c
@@ -19,6 +19,12 @@
 
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+
 /*
  * invpio2:  53 bits of 2/pi
  * pio2_1:   first  33 bit of pi/2
@@ -29,6 +35,7 @@
  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  */
 static const double
+toint   = 1.5/EPS,
 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 pio2_1  = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
@@ -41,8 +48,8 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 int __rem_pio2(double x, double *y)
 {
 	union {double f; uint64_t i;} u = {x};
-	double_t z,w,t,r;
-	double tx[3],ty[2],fn;
+	double_t z,w,t,r,fn;
+	double tx[3],ty[2];
 	uint32_t ix;
 	int sign, n, ex, ey, i;
 
@@ -111,8 +118,7 @@ int __rem_pio2(double x, double *y)
 	if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
 medium:
 		/* rint(x/(pi/2)), Assume round-to-nearest. */
-		fn = x*invpio2 + 0x1.8p52;
-		fn = fn - 0x1.8p52;
+		fn = x*invpio2 + toint - toint;
 		n = (int32_t)fn;
 		r = x - fn*pio2_1;
 		w = fn*pio2_1t;  /* 1st round, good to 85 bits */
diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c
index 838e1fca..f5163856 100644
--- a/src/math/__rem_pio2f.c
+++ b/src/math/__rem_pio2f.c
@@ -22,12 +22,19 @@
 
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+
 /*
  * invpio2:  53 bits of 2/pi
  * pio2_1:   first 25 bits of pi/2
  * pio2_1t:  pi/2 - pio2_1
  */
 static const double
+toint   = 1.5/EPS,
 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 pio2_1  = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
 pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
@@ -35,7 +42,8 @@ pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
 int __rem_pio2f(float x, double *y)
 {
 	union {float f; uint32_t i;} u = {x};
-	double tx[1],ty[1],fn;
+	double tx[1],ty[1];
+	double_t fn;
 	uint32_t ix;
 	int n, sign, e0;
 
@@ -43,8 +51,7 @@ int __rem_pio2f(float x, double *y)
 	/* 25+53 bit pi is good enough for medium size */
 	if (ix < 0x4dc90fdb) {  /* |x| ~< 2^28*(pi/2), medium size */
 		/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-		fn = x*invpio2 + 0x1.8p52;
-		fn = fn - 0x1.8p52;
+		fn = x*invpio2 + toint - toint;
 		n  = (int32_t)fn;
 		*y = x - fn*pio2_1 - fn*pio2_1t;
 		return n;
diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c
index 8b15b7b2..77255bd8 100644
--- a/src/math/__rem_pio2l.c
+++ b/src/math/__rem_pio2l.c
@@ -20,10 +20,11 @@
  * use __rem_pio2_large() for large x
  */
 
+static const long double toint = 1.5/LDBL_EPSILON;
+
 #if LDBL_MANT_DIG == 64
 /* u ~< 0x1p25*pi/2 */
 #define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.m>>48) < ((0x3fff + 25)<<16 | 0x921f>>1 | 0x8000))
-#define TOINT 0x1.8p63
 #define QUOBITS(x) ((uint32_t)(int32_t)x & 0x7fffffff)
 #define ROUND1 22
 #define ROUND2 61
@@ -50,7 +51,6 @@ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
 #elif LDBL_MANT_DIG == 113
 /* u ~< 0x1p45*pi/2 */
 #define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.top) < ((0x3fff + 45)<<16 | 0x921f))
-#define TOINT 0x1.8p112
 #define QUOBITS(x) ((uint32_t)(int64_t)x & 0x7fffffff)
 #define ROUND1 51
 #define ROUND2 119
@@ -77,7 +77,7 @@ int __rem_pio2l(long double x, long double *y)
 	ex = u.i.se & 0x7fff;
 	if (SMALL(u)) {
 		/* rint(x/(pi/2)), Assume round-to-nearest. */
-		fn = x*invpio2 + TOINT - TOINT;
+		fn = x*invpio2 + toint - toint;
 		n = QUOBITS(fn);
 		r = x-fn*pio2_1;
 		w = fn*pio2_1t;  /* 1st round good to 102/180 bits (ld80/ld128) */
diff --git a/src/math/ceil.c b/src/math/ceil.c
index 22dc224c..b13e6f2d 100644
--- a/src/math/ceil.c
+++ b/src/math/ceil.c
@@ -1,5 +1,12 @@
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+static const double_t toint = 1/EPS;
+
 double ceil(double x)
 {
 	union {double f; uint64_t i;} u = {x};
@@ -10,9 +17,9 @@ double ceil(double x)
 		return x;
 	/* y = int(x) - x, where int(x) is an integer neighbor of x */
 	if (u.i >> 63)
-		y = (double)(x - 0x1p52) + 0x1p52 - x;
+		y = x - toint + toint - x;
 	else
-		y = (double)(x + 0x1p52) - 0x1p52 - x;
+		y = x + toint - toint - x;
 	/* special case because of non-nearest rounding modes */
 	if (e <= 0x3ff-1) {
 		FORCE_EVAL(y);
diff --git a/src/math/ceill.c b/src/math/ceill.c
index a2cb0a7f..60a83020 100644
--- a/src/math/ceill.c
+++ b/src/math/ceill.c
@@ -6,11 +6,9 @@ long double ceill(long double x)
 	return ceil(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double ceill(long double x)
 {
 	union ldshape u = {x};
@@ -21,9 +19,9 @@ long double ceill(long double x)
 		return x;
 	/* y = int(x) - x, where int(x) is an integer neighbor of x */
 	if (u.i.se >> 15)
-		y = x - TOINT + TOINT - x;
+		y = x - toint + toint - x;
 	else
-		y = x + TOINT - TOINT - x;
+		y = x + toint - toint - x;
 	/* special case because of non-nearest rounding modes */
 	if (e <= 0x3fff-1) {
 		FORCE_EVAL(y);
diff --git a/src/math/floor.c b/src/math/floor.c
index ebc9fabe..14a31cd8 100644
--- a/src/math/floor.c
+++ b/src/math/floor.c
@@ -1,5 +1,12 @@
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+static const double_t toint = 1/EPS;
+
 double floor(double x)
 {
 	union {double f; uint64_t i;} u = {x};
@@ -10,9 +17,9 @@ double floor(double x)
 		return x;
 	/* y = int(x) - x, where int(x) is an integer neighbor of x */
 	if (u.i >> 63)
-		y = (double)(x - 0x1p52) + 0x1p52 - x;
+		y = x - toint + toint - x;
 	else
-		y = (double)(x + 0x1p52) - 0x1p52 - x;
+		y = x + toint - toint - x;
 	/* special case because of non-nearest rounding modes */
 	if (e <= 0x3ff-1) {
 		FORCE_EVAL(y);
diff --git a/src/math/floorl.c b/src/math/floorl.c
index 961f9e89..16aaec48 100644
--- a/src/math/floorl.c
+++ b/src/math/floorl.c
@@ -6,11 +6,9 @@ long double floorl(long double x)
 	return floor(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double floorl(long double x)
 {
 	union ldshape u = {x};
@@ -21,9 +19,9 @@ long double floorl(long double x)
 		return x;
 	/* y = int(x) - x, where int(x) is an integer neighbor of x */
 	if (u.i.se >> 15)
-		y = x - TOINT + TOINT - x;
+		y = x - toint + toint - x;
 	else
-		y = x + TOINT - TOINT - x;
+		y = x + toint - toint - x;
 	/* special case because of non-nearest rounding modes */
 	if (e <= 0x3fff-1) {
 		FORCE_EVAL(y);
diff --git a/src/math/modfl.c b/src/math/modfl.c
index 4b03a4be..a47b1924 100644
--- a/src/math/modfl.c
+++ b/src/math/modfl.c
@@ -11,11 +11,9 @@ long double modfl(long double x, long double *iptr)
 	return r;
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double modfl(long double x, long double *iptr)
 {
 	union ldshape u = {x};
@@ -40,7 +38,7 @@ long double modfl(long double x, long double *iptr)
 
 	/* raises spurious inexact */
 	absx = s ? -x : x;
-	y = absx + TOINT - TOINT - absx;
+	y = absx + toint - toint - absx;
 	if (y == 0) {
 		*iptr = x;
 		return s ? -0.0 : 0.0;
diff --git a/src/math/rintl.c b/src/math/rintl.c
index 26725073..374327db 100644
--- a/src/math/rintl.c
+++ b/src/math/rintl.c
@@ -6,11 +6,9 @@ long double rintl(long double x)
 	return rint(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double rintl(long double x)
 {
 	union ldshape u = {x};
@@ -21,9 +19,9 @@ long double rintl(long double x)
 	if (e >= 0x3fff+LDBL_MANT_DIG-1)
 		return x;
 	if (s)
-		y = x - TOINT + TOINT;
+		y = x - toint + toint;
 	else
-		y = x + TOINT - TOINT;
+		y = x + toint - toint;
 	if (y == 0)
 		return 0*x;
 	return y;
diff --git a/src/math/round.c b/src/math/round.c
index 4b38d1fd..130d58d2 100644
--- a/src/math/round.c
+++ b/src/math/round.c
@@ -1,5 +1,12 @@
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+static const double_t toint = 1/EPS;
+
 double round(double x)
 {
 	union {double f; uint64_t i;} u = {x};
@@ -12,10 +19,10 @@ double round(double x)
 		x = -x;
 	if (e < 0x3ff-1) {
 		/* raise inexact if x!=0 */
-		FORCE_EVAL(x + 0x1p52);
+		FORCE_EVAL(x + toint);
 		return 0*u.f;
 	}
-	y = (double)(x + 0x1p52) - 0x1p52 - x;
+	y = x + toint - toint - x;
 	if (y > 0.5)
 		y = y + x - 1;
 	else if (y <= -0.5)
diff --git a/src/math/roundf.c b/src/math/roundf.c
index c6b27797..e8210af5 100644
--- a/src/math/roundf.c
+++ b/src/math/roundf.c
@@ -1,5 +1,14 @@
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0
+#define EPS FLT_EPSILON
+#elif FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+static const float_t toint = 1/EPS;
+
 float roundf(float x)
 {
 	union {float f; uint32_t i;} u = {x};
@@ -11,10 +20,10 @@ float roundf(float x)
 	if (u.i >> 31)
 		x = -x;
 	if (e < 0x7f-1) {
-		FORCE_EVAL(x + 0x1p23f);
+		FORCE_EVAL(x + toint);
 		return 0*u.f;
 	}
-	y = (float)(x + 0x1p23f) - 0x1p23f - x;
+	y = x + toint - toint - x;
 	if (y > 0.5f)
 		y = y + x - 1;
 	else if (y <= -0.5f)
diff --git a/src/math/roundl.c b/src/math/roundl.c
index 8f3f28b3..f4ff6820 100644
--- a/src/math/roundl.c
+++ b/src/math/roundl.c
@@ -6,11 +6,9 @@ long double roundl(long double x)
 	return round(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double roundl(long double x)
 {
 	union ldshape u = {x};
@@ -22,10 +20,10 @@ long double roundl(long double x)
 	if (u.i.se >> 15)
 		x = -x;
 	if (e < 0x3fff-1) {
-		FORCE_EVAL(x + TOINT);
+		FORCE_EVAL(x + toint);
 		return 0*u.f;
 	}
-	y = x + TOINT - TOINT - x;
+	y = x + toint - toint - x;
 	if (y > 0.5)
 		y = y + x - 1;
 	else if (y <= -0.5)
diff --git a/src/math/truncl.c b/src/math/truncl.c
index 3eedb083..f07b1934 100644
--- a/src/math/truncl.c
+++ b/src/math/truncl.c
@@ -6,11 +6,9 @@ long double truncl(long double x)
 	return trunc(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double truncl(long double x)
 {
 	union ldshape u = {x};
@@ -27,7 +25,7 @@ long double truncl(long double x)
 	/* y = int(|x|) - |x|, where int(|x|) is an integer neighbor of |x| */
 	if (s)
 		x = -x;
-	y = x + TOINT - TOINT - x;
+	y = x + toint - toint - x;
 	if (y > 0)
 		y -= 1;
 	x += y;
generated by cgit-pink 1.4.1 (git 2.36.1) at 2024-10-31 06:06:57 +0000