From c3587effe27a3ac8c1406f064b7705963be9887a Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 19:14:32 +0100
Subject: minor fix in __tanl (get sign properly)

---
 src/math/__tanl.c | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

(limited to 'src/math')

diff --git a/src/math/__tanl.c b/src/math/__tanl.c
index e39e9df4..50ba2140 100644
--- a/src/math/__tanl.c
+++ b/src/math/__tanl.c
@@ -51,8 +51,7 @@ long double __tanl(long double x, long double y, int iy) {
 	int i;
 
 	iy = iy == 1 ? -1 : 1;        /* XXX recover original interface */
-	// FIXME: this is wrong, use copysign, signbit or union bithack
-	osign = x >= 0 ? 1.0 : -1.0;  /* XXX slow, probably wrong for -0 */
+	osign = copysignl(1.0, x);
 	if (fabsl(x) >= 0.67434) {
 		if (x < 0) {
 			x = -x;
-- 
cgit v1.2.1


From 2e8c8fbe7d65ba0026cb084dc8570d94cbc908ff Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 19:26:31 +0100
Subject: don't inline __rem_pio2l so the code size is smaller

---
 src/math/__rem_pio2l.c | 131 +++++++++++++++++++++++++++++++++++++++++++++++++
 src/math/__rem_pio2l.h | 131 -------------------------------------------------
 src/math/cosl.c        |   2 -
 src/math/sincosl.c     |   2 -
 src/math/sinl.c        |   2 -
 src/math/tanl.c        |   2 -
 6 files changed, 131 insertions(+), 139 deletions(-)
 create mode 100644 src/math/__rem_pio2l.c
 delete mode 100644 src/math/__rem_pio2l.h

(limited to 'src/math')

diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c
new file mode 100644
index 00000000..10af404c
--- /dev/null
+++ b/src/math/__rem_pio2l.c
@@ -0,0 +1,131 @@
+/* origin: FreeBSD /usr/src/lib/msun/ld80/e_rem_pio2.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+#include "libm.h"
+#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+/* ld80 version of __rem_pio2(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __rem_pio2_large() for large x
+ */
+
+#define BIAS    (LDBL_MAX_EXP - 1)
+
+/*
+ * invpio2:  64 bits of 2/pi
+ * pio2_1:   first  39 bits of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 39 bits of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  39 bits of pi/2
+ * 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 */
+pio2_3 =  6.36831716351370313614e-25; /*  0x18a2e037074000.0p-133 */
+
+static const long double
+invpio2 =  6.36619772367581343076e-01L, /*  0xa2f9836e4e44152a.0p-64 */
+pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
+pio2_2t =  6.36831716351095013979e-25L, /*  0xc51701b839a25205.0p-144 */
+pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
+
+int __rem_pio2l(long double x, long double *y)
+{
+	union IEEEl2bits u,u1;
+	long double z,w,t,r,fn;
+	double tx[3],ty[2];
+	int e0,ex,i,j,nx,n;
+	int16_t expsign;
+
+	u.e = x;
+	expsign = u.xbits.expsign;
+	ex = expsign & 0x7fff;
+	if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) {
+		union IEEEl2bits u2;
+		int ex1;
+
+		/* |x| ~< 2^25*(pi/2), medium size */
+		/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+		fn = x*invpio2 + 0x1.8p63;
+		fn = fn - 0x1.8p63;
+// FIXME
+//#ifdef HAVE_EFFICIENT_IRINT
+//		n = irint(fn);
+//#else
+		n = fn;
+//#endif
+		r = x-fn*pio2_1;
+		w = fn*pio2_1t;    /* 1st round good to 102 bit */
+		j = ex;
+		y[0] = r-w;
+		u2.e = y[0];
+		ex1 = u2.xbits.expsign & 0x7fff;
+		i = j-ex1;
+		if (i > 22) {  /* 2nd iteration needed, good to 141 */
+			t = r;
+			w = fn*pio2_2;
+			r = t-w;
+			w = fn*pio2_2t-((t-r)-w);
+			y[0] = r-w;
+			u2.e = y[0];
+			ex1 = u2.xbits.expsign & 0x7fff;
+			i = j-ex1;
+			if (i > 61) {  /* 3rd iteration need, 180 bits acc */
+				t = r; /* will cover all possible cases */
+				w = fn*pio2_3;
+				r = t-w;
+				w = fn*pio2_3t-((t-r)-w);
+				y[0] = r-w;
+			}
+		}
+		y[1] = (r - y[0]) - w;
+		return n;
+	}
+	/*
+	 * all other (large) arguments
+	 */
+	if (ex == 0x7fff) {                /* x is inf or NaN */
+		y[0] = y[1] = x - x;
+		return 0;
+	}
+	/* set z = scalbn(|x|,ilogb(x)-23) */
+	u1.e = x;
+	e0 = ex - BIAS - 23;            /* e0 = ilogb(|x|)-23; */
+	u1.xbits.expsign = ex - e0;
+	z = u1.e;
+	for (i=0; i<2; i++) {
+		tx[i] = (double)(int32_t)z;
+		z     = (z-tx[i])*two24;
+	}
+	tx[2] = z;
+	nx = 3;
+	while (tx[nx-1] == zero)
+		nx--;     /* skip zero term */
+	n = __rem_pio2_large(tx,ty,e0,nx,2);
+	r = (long double)ty[0] + ty[1];
+	w = ty[1] - (r - ty[0]);
+	if (expsign < 0) {
+		y[0] = -r;
+		y[1] = -w;
+		return -n;
+	}
+	y[0] = r;
+	y[1] = w;
+	return n;
+}
+#endif
diff --git a/src/math/__rem_pio2l.h b/src/math/__rem_pio2l.h
deleted file mode 100644
index 11123c3f..00000000
--- a/src/math/__rem_pio2l.h
+++ /dev/null
@@ -1,131 +0,0 @@
-/* origin: FreeBSD /usr/src/lib/msun/ld80/e_rem_pio2.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
- *
- * 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.
- * ====================================================
- *
- * Optimized by Bruce D. Evans.
- */
-#include "libm.h"
-#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-/* ld80 version of __rem_pio2(x,y)
- *
- * return the remainder of x rem pi/2 in y[0]+y[1]
- * use __rem_pio2_large() for large x
- */
-
-#define BIAS    (LDBL_MAX_EXP - 1)
-
-/*
- * invpio2:  64 bits of 2/pi
- * pio2_1:   first  39 bits of pi/2
- * pio2_1t:  pi/2 - pio2_1
- * pio2_2:   second 39 bits of pi/2
- * pio2_2t:  pi/2 - (pio2_1+pio2_2)
- * pio2_3:   third  39 bits of pi/2
- * 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 */
-pio2_3 =  6.36831716351370313614e-25; /*  0x18a2e037074000.0p-133 */
-
-static const long double
-invpio2 =  6.36619772367581343076e-01L, /*  0xa2f9836e4e44152a.0p-64 */
-pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
-pio2_2t =  6.36831716351095013979e-25L, /*  0xc51701b839a25205.0p-144 */
-pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
-
-static inline int __rem_pio2l(long double x, long double *y)
-{
-	union IEEEl2bits u,u1;
-	long double z,w,t,r,fn;
-	double tx[3],ty[2];
-	int e0,ex,i,j,nx,n;
-	int16_t expsign;
-
-	u.e = x;
-	expsign = u.xbits.expsign;
-	ex = expsign & 0x7fff;
-	if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) {
-		union IEEEl2bits u2;
-		int ex1;
-
-		/* |x| ~< 2^25*(pi/2), medium size */
-		/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-		fn = x*invpio2 + 0x1.8p63;
-		fn = fn - 0x1.8p63;
-// FIXME
-//#ifdef HAVE_EFFICIENT_IRINT
-//		n = irint(fn);
-//#else
-		n = fn;
-//#endif
-		r = x-fn*pio2_1;
-		w = fn*pio2_1t;    /* 1st round good to 102 bit */
-		j = ex;
-		y[0] = r-w;
-		u2.e = y[0];
-		ex1 = u2.xbits.expsign & 0x7fff;
-		i = j-ex1;
-		if (i > 22) {  /* 2nd iteration needed, good to 141 */
-			t = r;
-			w = fn*pio2_2;
-			r = t-w;
-			w = fn*pio2_2t-((t-r)-w);
-			y[0] = r-w;
-			u2.e = y[0];
-			ex1 = u2.xbits.expsign & 0x7fff;
-			i = j-ex1;
-			if (i > 61) {  /* 3rd iteration need, 180 bits acc */
-				t = r; /* will cover all possible cases */
-				w = fn*pio2_3;
-				r = t-w;
-				w = fn*pio2_3t-((t-r)-w);
-				y[0] = r-w;
-			}
-		}
-		y[1] = (r - y[0]) - w;
-		return n;
-	}
-	/*
-	 * all other (large) arguments
-	 */
-	if (ex == 0x7fff) {                /* x is inf or NaN */
-		y[0] = y[1] = x - x;
-		return 0;
-	}
-	/* set z = scalbn(|x|,ilogb(x)-23) */
-	u1.e = x;
-	e0 = ex - BIAS - 23;            /* e0 = ilogb(|x|)-23; */
-	u1.xbits.expsign = ex - e0;
-	z = u1.e;
-	for (i=0; i<2; i++) {
-		tx[i] = (double)(int32_t)z;
-		z     = (z-tx[i])*two24;
-	}
-	tx[2] = z;
-	nx = 3;
-	while (tx[nx-1] == zero)
-		nx--;     /* skip zero term */
-	n = __rem_pio2_large(tx,ty,e0,nx,2);
-	r = (long double)ty[0] + ty[1];
-	w = ty[1] - (r - ty[0]);
-	if (expsign < 0) {
-		y[0] = -r;
-		y[1] = -w;
-		return -n;
-	}
-	y[0] = r;
-	y[1] = w;
-	return n;
-}
-#endif
diff --git a/src/math/cosl.c b/src/math/cosl.c
index 2c650cdc..25d9005a 100644
--- a/src/math/cosl.c
+++ b/src/math/cosl.c
@@ -36,8 +36,6 @@ long double cosl(long double x) {
 	return cos(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#include "__rem_pio2l.h"
-
 long double cosl(long double x)
 {
 	union IEEEl2bits z;
diff --git a/src/math/sincosl.c b/src/math/sincosl.c
index 378dc979..e14129a2 100644
--- a/src/math/sincosl.c
+++ b/src/math/sincosl.c
@@ -9,8 +9,6 @@ void sincosl(long double x, long double *sin, long double *cos)
 	*cos = c;
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#include "__rem_pio2l.h"
-
 void sincosl(long double x, long double *sin, long double *cos)
 {
 	union IEEEl2bits u;
diff --git a/src/math/sinl.c b/src/math/sinl.c
index 0b1aeb75..7e0b44f4 100644
--- a/src/math/sinl.c
+++ b/src/math/sinl.c
@@ -34,8 +34,6 @@ long double sinl(long double x)
 	return sin(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#include "__rem_pio2l.h"
-
 long double sinl(long double x)
 {
 	union IEEEl2bits z;
diff --git a/src/math/tanl.c b/src/math/tanl.c
index 462ead91..0194eaf7 100644
--- a/src/math/tanl.c
+++ b/src/math/tanl.c
@@ -38,8 +38,6 @@ long double tanl(long double x)
 	return tan(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#include "__rem_pio2l.h"
-
 long double tanl(long double x)
 {
 	union IEEEl2bits z;
-- 
cgit v1.2.1


From 01fdfd491b5d83b72099fbae14c4a71ed8e0b945 Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 22:49:03 +0100
Subject: fix long double const workaround in cbrtl

---
 src/math/cbrtl.c | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

(limited to 'src/math')

diff --git a/src/math/cbrtl.c b/src/math/cbrtl.c
index d138b9f2..5297d68f 100644
--- a/src/math/cbrtl.c
+++ b/src/math/cbrtl.c
@@ -118,11 +118,7 @@ long double cbrtl(long double x)
 	 * Round it away from zero to 32 bits (32 so that t*t is exact, and
 	 * away from zero for technical reasons).
 	 */
-	volatile double vd2 = 0x1.0p32;
-	volatile double vd1 = 0x1.0p-31;
-	#define vd ((long double)vd2 + vd1)
-
-	t = dt + vd - 0x1.0p32;
+	t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32;
 #elif LDBL_MANT_DIG == 113
 	/*
 	 * Round dt away from zero to 47 bits.  Since we don't trust the 47,
-- 
cgit v1.2.1


From 2786c7d21611b9fa3b2fe356542cf213e7dd0ba4 Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 22:57:58 +0100
Subject: use scalbn or *2.0 instead of ldexp, fix fmal

Some code assumed ldexp(x, 1) is faster than 2.0*x,
but ldexp is a wrapper around scalbn which uses
multiplications inside, so this optimization is
wrong.

This commit also fixes fmal which accidentally
used ldexp instead of ldexpl loosing precision.

There are various additional changes from the
work-in-progress const cleanups.
---
 src/math/expl.c   |   8 ++---
 src/math/expm1l.c |   6 ++--
 src/math/fma.c    |  12 +++----
 src/math/fmal.c   |  12 +++----
 src/math/log10l.c |  22 ++++++------
 src/math/log2l.c  |  20 +++++------
 src/math/logl.c   |  20 +++++------
 src/math/powl.c   | 103 +++++++++++++++++++++++++++---------------------------
 8 files changed, 102 insertions(+), 101 deletions(-)

(limited to 'src/math')

diff --git a/src/math/expl.c b/src/math/expl.c
index 9507fd2e..b289ffec 100644
--- a/src/math/expl.c
+++ b/src/math/expl.c
@@ -102,13 +102,13 @@ long double expl(long double x)
 	if (x > MAXLOGL)
 		return INFINITY;
 	if (x < MINLOGL)
-		return 0.0L;
+		return 0.0;
 
 	/* Express e**x = e**g 2**n
 	 *   = e**g e**(n loge(2))
 	 *   = e**(g + n loge(2))
 	 */
-	px = floorl(LOG2EL * x + 0.5L); /* floor() truncates toward -infinity. */
+	px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */
 	n = px;
 	x -= px * C1;
 	x -= px * C2;
@@ -120,8 +120,8 @@ long double expl(long double x)
 	xx = x * x;
 	px = x * __polevll(xx, P, 2);
 	x =  px/(__polevll(xx, Q, 3) - px);
-	x = 1.0L + ldexpl(x, 1);
-	x = ldexpl(x, n);
+	x = 1.0 + 2.0 * x;
+	x = scalbnl(x, n);
 	return x;
 }
 #endif
diff --git a/src/math/expm1l.c b/src/math/expm1l.c
index 2f94dfa2..ca0d141f 100644
--- a/src/math/expm1l.c
+++ b/src/math/expm1l.c
@@ -97,11 +97,11 @@ long double expm1l(long double x)
 		return x;
 	/* Minimum value.*/
 	if (x < minarg)
-		return -1.0L;
+		return -1.0;
 
 	xx = C1 + C2;
 	/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
-	px = floorl (0.5 + x / xx);
+	px = floorl(0.5 + x / xx);
 	k = px;
 	/* remainder times ln 2 */
 	x -= px * C1;
@@ -116,7 +116,7 @@ long double expm1l(long double x)
 	/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 	 We have qx = exp(remainder ln 2) - 1, so
 	 exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.  */
-	px = ldexpl(1.0L, k);
+	px = scalbnl(1.0, k);
 	x = px * qx + (px - 1.0);
 	return x;
 }
diff --git a/src/math/fma.c b/src/math/fma.c
index 87d450c7..5fb95406 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -247,7 +247,7 @@ static inline double add_and_denormalize(double a, double b, int scale)
 			INSERT_WORD64(sum.hi, hibits);
 		}
 	}
-	return (ldexp(sum.hi, scale));
+	return scalbn(sum.hi, scale);
 }
 
 /*
@@ -364,7 +364,7 @@ double fma(double x, double y, double z)
 		}
 	}
 	if (spread <= DBL_MANT_DIG * 2)
-		zs = ldexp(zs, -spread);
+		zs = scalbn(zs, -spread);
 	else
 		zs = copysign(DBL_MIN, zs);
 
@@ -390,7 +390,7 @@ double fma(double x, double y, double z)
 		 */
 		fesetround(oround);
 		volatile double vzs = zs; /* XXX gcc CSE bug workaround */
-		return (xy.hi + vzs + ldexp(xy.lo, spread));
+		return xy.hi + vzs + scalbn(xy.lo, spread);
 	}
 
 	if (oround != FE_TONEAREST) {
@@ -400,13 +400,13 @@ double fma(double x, double y, double z)
 		 */
 		fesetround(oround);
 		adj = r.lo + xy.lo;
-		return (ldexp(r.hi + adj, spread));
+		return scalbn(r.hi + adj, spread);
 	}
 
 	adj = add_adjusted(r.lo, xy.lo);
 	if (spread + ilogb(r.hi) > -1023)
-		return (ldexp(r.hi + adj, spread));
+		return scalbn(r.hi + adj, spread);
 	else
-		return (add_and_denormalize(r.hi, adj, spread));
+		return add_and_denormalize(r.hi, adj, spread);
 }
 #endif
diff --git a/src/math/fmal.c b/src/math/fmal.c
index cbaf46eb..be64f145 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -115,7 +115,7 @@ static inline long double add_and_denormalize(long double a, long double b, int
 		if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
 			sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
 	}
-	return (ldexp(sum.hi, scale));
+	return scalbnl(sum.hi, scale);
 }
 
 /*
@@ -228,7 +228,7 @@ long double fmal(long double x, long double y, long double z)
 		}
 	}
 	if (spread <= LDBL_MANT_DIG * 2)
-		zs = ldexpl(zs, -spread);
+		zs = scalbnl(zs, -spread);
 	else
 		zs = copysignl(LDBL_MIN, zs);
 
@@ -254,7 +254,7 @@ long double fmal(long double x, long double y, long double z)
 		 */
 		fesetround(oround);
 		volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
-		return (xy.hi + vzs + ldexpl(xy.lo, spread));
+		return xy.hi + vzs + scalbnl(xy.lo, spread);
 	}
 
 	if (oround != FE_TONEAREST) {
@@ -264,13 +264,13 @@ long double fmal(long double x, long double y, long double z)
 		 */
 		fesetround(oround);
 		adj = r.lo + xy.lo;
-		return (ldexpl(r.hi + adj, spread));
+		return scalbnl(r.hi + adj, spread);
 	}
 
 	adj = add_adjusted(r.lo, xy.lo);
 	if (spread + ilogbl(r.hi) > -16383)
-		return (ldexpl(r.hi + adj, spread));
+		return scalbnl(r.hi + adj, spread);
 	else
-		return (add_and_denormalize(r.hi, adj, spread));
+		return add_and_denormalize(r.hi, adj, spread);
 }
 #endif
diff --git a/src/math/log10l.c b/src/math/log10l.c
index b954cc77..f0eeeafb 100644
--- a/src/math/log10l.c
+++ b/src/math/log10l.c
@@ -123,9 +123,9 @@ long double log10l(long double x)
 
 	if (isnan(x))
 		return x;
-	if(x <= 0.0L) {
-		if(x == 0.0L)
-			return -1.0L / (x - x);
+	if(x <= 0.0) {
+		if(x == 0.0)
+			return -1.0 / (x - x);
 		return (x - x) / (x - x);
 	}
 	if (x == INFINITY)
@@ -142,12 +142,12 @@ long double log10l(long double x)
 	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)   */
-			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;
 		}
 		x = z / y;
 		z = x*x;
@@ -158,13 +158,13 @@ long double log10l(long double x)
 	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 	if (x < SQRTH) {
 		e -= 1;
-		x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+		x = 2.0*x - 1.0;
 	} else {
-		x = x - 1.0L;
+		x = x - 1.0;
 	}
 	z = x*x;
 	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
-	y = y - ldexpl(z, -1);   /* -0.5x^2 + ... */
+	y = y - 0.5*z;
 
 done:
 	/* Multiply log of fraction by log10(e)
diff --git a/src/math/log2l.c b/src/math/log2l.c
index 4339c033..8ebce9c4 100644
--- a/src/math/log2l.c
+++ b/src/math/log2l.c
@@ -121,8 +121,8 @@ long double log2l(long double x)
 		return x;
 	if (x == INFINITY)
 		return x;
-	if (x <= 0.0L) {
-		if (x == 0.0L)
+	if (x <= 0.0) {
+		if (x == 0.0)
 			return -INFINITY;
 		return NAN;
 	}
@@ -139,12 +139,12 @@ long double log2l(long double x)
 	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)   */
-			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;
 		}
 		x = z / y;
 		z = x*x;
@@ -155,13 +155,13 @@ long double log2l(long double x)
 	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 	if (x < SQRTH) {
 		e -= 1;
-		x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+		x = 2.0*x - 1.0;
 	} else {
-		x = x - 1.0L;
+		x = x - 1.0;
 	}
 	z = x*x;
 	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
-	y = y - ldexpl(z, -1);   /* -0.5x^2 + ... */
+	y = y - 0.5*z;
 
 done:
 	/* Multiply log of fraction by log2(e)
diff --git a/src/math/logl.c b/src/math/logl.c
index ee7ca64a..ffd81038 100644
--- a/src/math/logl.c
+++ b/src/math/logl.c
@@ -119,8 +119,8 @@ long double logl(long double x)
 		return x;
 	if (x == INFINITY)
 		return x;
-	if (x <= 0.0L) {
-		if (x == 0.0L)
+	if (x <= 0.0) {
+		if (x == 0.0)
 			return -INFINITY;
 		return NAN;
 	}
@@ -137,12 +137,12 @@ long double logl(long double x)
 	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)   */
-			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;
 		}
 		x = z / y;
 		z = x*x;
@@ -156,14 +156,14 @@ long double logl(long double x)
 	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 	if (x < SQRTH) {
 		e -= 1;
-		x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+		x = 2.0*x - 1.0;
 	} else {
-		x = x - 1.0L;
+		x = x - 1.0;
 	}
 	z = x*x;
 	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
 	y = y + e * C2;
-	z = y - ldexpl(z, -1);   /*  y - 0.5 * z  */
+	z = y - 0.5*z;
 	/* Note, the sum of above terms does not exceed x/4,
 	 * so it contributes at most about 1/4 lsb to the error.
 	 */
diff --git a/src/math/powl.c b/src/math/powl.c
index a6ee275f..a1d2f076 100644
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -203,44 +203,44 @@ long double powl(long double x, long double y)
 	volatile long double z=0;
 	long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
 
-	if (y == 0.0L)
-		return 1.0L;
+	if (y == 0.0)
+		return 1.0;
 	if (isnan(x))
 		return x;
 	if (isnan(y))
 		return y;
-	if (y == 1.0L)
+	if (y == 1.0)
 		return x;
 
 	// FIXME: this is wrong, see pow special cases in c99 F.9.4.4
-	if (!isfinite(y) && (x == -1.0L || x == 1.0L) )
+	if (!isfinite(y) && (x == -1.0 || x == 1.0) )
 		return y - y;   /* +-1**inf is NaN */
-	if (x == 1.0L)
-		return 1.0L;
+	if (x == 1.0)
+		return 1.0;
 	if (y >= LDBL_MAX) {
-		if (x > 1.0L)
+		if (x > 1.0)
 			return INFINITY;
-		if (x > 0.0L && x < 1.0L)
-			return 0.0L;
-		if (x < -1.0L)
+		if (x > 0.0 && x < 1.0)
+			return 0.0;
+		if (x < -1.0)
 			return INFINITY;
-		if (x > -1.0L && x < 0.0L)
-			return 0.0L;
+		if (x > -1.0 && x < 0.0)
+			return 0.0;
 	}
 	if (y <= -LDBL_MAX) {
-		if (x > 1.0L)
-			return 0.0L;
-		if (x > 0.0L && x < 1.0L)
+		if (x > 1.0)
+			return 0.0;
+		if (x > 0.0 && x < 1.0)
 			return INFINITY;
-		if (x < -1.0L)
-			return 0.0L;
-		if (x > -1.0L && x < 0.0L)
+		if (x < -1.0)
+			return 0.0;
+		if (x > -1.0 && x < 0.0)
 			return INFINITY;
 	}
 	if (x >= LDBL_MAX) {
-		if (y > 0.0L)
+		if (y > 0.0)
 			return INFINITY;
-		return 0.0L;
+		return 0.0;
 	}
 
 	w = floorl(y);
@@ -253,29 +253,29 @@ long double powl(long double x, long double y)
 	yoddint = 0;
 	if (iyflg) {
 		ya = fabsl(y);
-		ya = floorl(0.5L * ya);
-		yb = 0.5L * fabsl(w);
+		ya = floorl(0.5 * ya);
+		yb = 0.5 * fabsl(w);
 		if( ya != yb )
 			yoddint = 1;
 	}
 
 	if (x <= -LDBL_MAX) {
-		if (y > 0.0L) {
+		if (y > 0.0) {
 			if (yoddint)
 				return -INFINITY;
 			return INFINITY;
 		}
-		if (y < 0.0L) {
+		if (y < 0.0) {
 			if (yoddint)
-				return -0.0L;
+				return -0.0;
 			return 0.0;
 		}
 	}
 
 
 	nflg = 0;       /* flag = 1 if x<0 raised to integer power */
-	if (x <= 0.0L) {
-		if (x == 0.0L) {
+	if (x <= 0.0) {
+		if (x == 0.0) {
 			if (y < 0.0) {
 				if (signbit(x) && yoddint)
 					return -INFINITY;
@@ -283,12 +283,12 @@ long double powl(long double x, long double y)
 			}
 			if (y > 0.0) {
 				if (signbit(x) && yoddint)
-					return -0.0L;
+					return -0.0;
 				return 0.0;
 			}
-			if (y == 0.0L)
-				return 1.0L;  /*   0**0   */
-			return 0.0L;  /*   0**y   */
+			if (y == 0.0)
+				return 1.0;  /*   0**0   */
+			return 0.0;  /*   0**y   */
 		}
 		if (iyflg == 0)
 			return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
@@ -343,7 +343,7 @@ long double powl(long double x, long double y)
 	 */
 	z = x*x;
 	w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
-	w = w - ldexpl(z, -1);  /*  w - 0.5 * z  */
+	w = w - 0.5*z;
 
 	/* Convert to base 2 logarithm:
 	 * multiply by log2(e) = 1 + LOG2EA
@@ -355,7 +355,8 @@ long double powl(long double x, long double y)
 
 	/* Compute exponent term of the base 2 logarithm. */
 	w = -i;
-	w = ldexpl(w, -LNXT); /* divide by NXT */
+	// TODO: use w * 0x1p-5;
+	w = scalbnl(w, -LNXT); /* divide by NXT */
 	w += e;
 	/* Now base 2 log of x is w + z. */
 
@@ -380,7 +381,7 @@ long double powl(long double x, long double y)
 
 	H = Fb + Gb;
 	Ha = reducl(H);
-	w = ldexpl( Ga+Ha, LNXT );
+	w = scalbnl( Ga+Ha, LNXT );
 
 	/* Test the power of 2 for overflow */
 	if (w > MEXP)
@@ -391,9 +392,9 @@ long double powl(long double x, long double y)
 	e = w;
 	Hb = H - Ha;
 
-	if (Hb > 0.0L) {
+	if (Hb > 0.0) {
 		e += 1;
-		Hb -= 1.0L/NXT;  /*0.0625L;*/
+		Hb -= 1.0/NXT;  /*0.0625L;*/
 	}
 
 	/* Now the product y * log2(x)  =  Hb + e/NXT.
@@ -415,16 +416,16 @@ long double powl(long double x, long double y)
 	w = douba(e);
 	z = w * z;  /*  2**-e * ( 1 + (2**Hb-1) )  */
 	z = z + w;
-	z = ldexpl(z, i);  /* multiply by integer power of 2 */
+	z = scalbnl(z, i);  /* multiply by integer power of 2 */
 
 	if (nflg) {
 		/* For negative x,
 		 * find out if the integer exponent
 		 * is odd or even.
 		 */
-		w = ldexpl(y, -1);
+		w = 0.5*y;
 		w = floorl(w);
-		w = ldexpl(w, 1);
+		w = 2.0*w;
 		if (w != y)
 			z = -z;  /* odd exponent */
 	}
@@ -438,9 +439,9 @@ static long double reducl(long double x)
 {
 	long double t;
 
-	t = ldexpl(x, LNXT);
+	t = scalbnl(x, LNXT);
 	t = floorl(t);
-	t = ldexpl(t, -LNXT);
+	t = scalbnl(t, -LNXT);
 	return t;
 }
 
@@ -483,18 +484,18 @@ static long double powil(long double x, int nn)
 	long double s;
 	int n, e, sign, asign, lx;
 
-	if (x == 0.0L) {
+	if (x == 0.0) {
 		if (nn == 0)
-			return 1.0L;
+			return 1.0;
 		else if (nn < 0)
 			return LDBL_MAX;
-		return 0.0L;
+		return 0.0;
 	}
 
 	if (nn == 0)
-		return 1.0L;
+		return 1.0;
 
-	if (x < 0.0L) {
+	if (x < 0.0) {
 		asign = -1;
 		x = -x;
 	} else
@@ -516,7 +517,7 @@ static long double powil(long double x, int nn)
 	e = (lx - 1)*n;
 	if ((e == 0) || (e > 64) || (e < -64)) {
 		s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);
-		s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
+		s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L;
 	} else {
 		s = LOGE2L * e;
 	}
@@ -530,8 +531,8 @@ static long double powil(long double x, int nn)
 	 * since roundoff error in 1.0/x will be amplified.
 	 * The precise demarcation should be the gradual underflow threshold.
 	 */
-	if (s < -MAXLOGL+2.0L) {
-		x = 1.0L/x;
+	if (s < -MAXLOGL+2.0) {
+		x = 1.0/x;
 		sign = -sign;
 	}
 
@@ -539,7 +540,7 @@ static long double powil(long double x, int nn)
 	if (n & 1)
 		y = x;
 	else {
-		y = 1.0L;
+		y = 1.0;
 		asign = 0;
 	}
 
@@ -555,7 +556,7 @@ static long double powil(long double x, int nn)
 	if (asign)
 		y = -y;  /* odd power of negative number */
 	if (sign < 0)
-		y = 1.0L/y;
+		y = 1.0/y;
 	return y;
 }
 
-- 
cgit v1.2.1


From 75483499dad38b97f5dabb710635e6a8bbbb1c84 Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 23:27:45 +0100
Subject: new modff.c code, fix nan handling in modfl

---
 src/math/modff.c | 66 +++++++++++++++++++++-----------------------------------
 src/math/modfl.c |  4 ++--
 2 files changed, 26 insertions(+), 44 deletions(-)

(limited to 'src/math')

diff --git a/src/math/modff.c b/src/math/modff.c
index d535314c..bf6e4ced 100644
--- a/src/math/modff.c
+++ b/src/math/modff.c
@@ -1,51 +1,33 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_modff.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 one = 1.0;
-
 float modff(float x, float *iptr)
 {
-	int32_t i0,j0;
-	uint32_t i;
+	uint32_t u, mask;
+	int e;
 
-	GET_FLOAT_WORD(i0, x);
-	j0 = ((i0>>23) & 0xff) - 0x7f;  /* exponent of x */
-	if (j0 < 23) {  /* integer part in x */
-		if (j0 < 0) {  /* |x| < 1 */
-			SET_FLOAT_WORD(*iptr, i0 & 0x80000000);  /* *iptr = +-0 */
-			return x;
-		}
-		i = 0x007fffff >> j0;
-		if ((i0&i) == 0) {  /* x is integral */
-			uint32_t ix;
-			*iptr = x;
-			GET_FLOAT_WORD(ix, x);
-			SET_FLOAT_WORD(x, ix & 0x80000000);  /* return +-0 */
-			return x;
-		}
-		SET_FLOAT_WORD(*iptr, i0&~i);
-		return x - *iptr;
-	} else {        /* no fraction part */
-		uint32_t ix;
-		*iptr = x*one;
-		if (x != x)  /* NaN */
+	GET_FLOAT_WORD(u, x);
+	e = (int)(u>>23 & 0xff) - 0x7f;
+
+	/* no fractional part */
+	if (e >= 23) {
+		*iptr = x;
+		if (e == 0x80 && u<<9 != 0) /* nan */
 			return x;
-		GET_FLOAT_WORD(ix, x);
-		SET_FLOAT_WORD(x, ix & 0x80000000);  /* return +-0 */
+		SET_FLOAT_WORD(x, u & 0x80000000);
+		return x;
+	}
+	/* no integral part */
+	if (e < 0) {
+		SET_FLOAT_WORD(*iptr, u & 0x80000000);
+		return x;
+	}
+
+	mask = 0x007fffff>>e;
+	if ((u & mask) == 0) {
+		*iptr = x;
+		SET_FLOAT_WORD(x, u & 0x80000000);
 		return x;
 	}
+	SET_FLOAT_WORD(*iptr, u & ~mask);
+	return x - *iptr;
 }
diff --git a/src/math/modfl.c b/src/math/modfl.c
index 2ca67b11..6520a1c2 100644
--- a/src/math/modfl.c
+++ b/src/math/modfl.c
@@ -54,7 +54,7 @@ long double modfl(long double x, long double *iptr)
 /* The number of fraction bits in manh, not counting the integer bit */
 #define HIBITS  (LDBL_MANT_DIG - LDBL_MANL_SIZE)
 
-static const long double zero[] = { 0.0L, -0.0L };
+static const long double zero[] = { 0.0, -0.0 };
 
 long double modfl(long double x, long double *iptr)
 {
@@ -81,7 +81,7 @@ long double modfl(long double x, long double *iptr)
 		return x - ux.e;
 	} else if (e >= LDBL_MANT_DIG - 1) {    /* x has no fraction part. */
 		*iptr = x;
-		if (x != x)                     /* Handle NaNs. */
+		if (e == LDBL_MAX_EXP && (ux.bits.manh|ux.bits.manl)) /* nan */
 			return x;
 		return zero[ux.bits.sign];
 	} else {                                /* Fraction part is in manl. */
-- 
cgit v1.2.1


From 4caa17b2a17d136efedfb63fceef832401063d70 Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 23:30:45 +0100
Subject: don't try to create non-standard denormalization signal

Underflow exception is only raised when the result is
invalid, but fmod is always exact. x87 has a denormalization
exception, but that's nonstandard. And the superflous *1.0
will be optimized away by any compiler that does not honor
signaling nans.
---
 src/math/fmod.c  | 3 +--
 src/math/fmodf.c | 3 +--
 src/math/fmodl.c | 5 ++---
 3 files changed, 4 insertions(+), 7 deletions(-)

(limited to 'src/math')

diff --git a/src/math/fmod.c b/src/math/fmod.c
index 6856844e..84a1b4ac 100644
--- a/src/math/fmod.c
+++ b/src/math/fmod.c
@@ -17,7 +17,7 @@
 
 #include "libm.h"
 
-static const double one = 1.0, Zero[] = {0.0, -0.0,};
+static const double Zero[] = {0.0, -0.0,};
 
 double fmod(double x, double y)
 {
@@ -140,7 +140,6 @@ double fmod(double x, double y)
 			lx = hx>>(n-32); hx = sx;
 		}
 		INSERT_WORDS(x, hx|sx, lx);
-		x *= one;  /* create necessary signal */
 	}
 	return x;  /* exact output */
 }
diff --git a/src/math/fmodf.c b/src/math/fmodf.c
index 4b50a3d3..837e9371 100644
--- a/src/math/fmodf.c
+++ b/src/math/fmodf.c
@@ -20,7 +20,7 @@
 
 #include "libm.h"
 
-static const float one = 1.0, Zero[] = {0.0, -0.0,};
+static const float Zero[] = {0.0, -0.0,};
 
 float fmodf(float x, float y)
 {
@@ -99,7 +99,6 @@ float fmodf(float x, float y)
 		n = -126 - iy;
 		hx >>= n;
 		SET_FLOAT_WORD(x, hx|sx);
-		x *= one;          /* create necessary signal */
 	}
 	return x;  /* exact output */
 }
diff --git a/src/math/fmodl.c b/src/math/fmodl.c
index 2e3eec1f..b930c49d 100644
--- a/src/math/fmodl.c
+++ b/src/math/fmodl.c
@@ -48,7 +48,7 @@ typedef uint32_t manh_t;
 
 #define MANL_SHIFT      (LDBL_MANL_SIZE - 1)
 
-static const long double one = 1.0, Zero[] = {0.0, -0.0,};
+static const long double Zero[] = {0.0, -0.0,};
 
 /*
  * fmodl(x,y)
@@ -153,7 +153,6 @@ long double fmodl(long double x, long double y)
 	} else {
 		ux.bits.exp = iy + BIAS;
 	}
-	x = ux.e * one;   /* create necessary signal */
-	return x;         /* exact output */
+	return ux.e;       /* exact output */
 }
 #endif
-- 
cgit v1.2.1


From b03255af77776703c8d48819e824d09f6f54ba86 Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 23:39:47 +0100
Subject: fix remainder*.c: remove useless long double cast

---
 src/math/remainder.c  | 12 ++++--------
 src/math/remainderf.c |  9 ++-------
 2 files changed, 6 insertions(+), 15 deletions(-)

(limited to 'src/math')

diff --git a/src/math/remainder.c b/src/math/remainder.c
index db176c88..2dede3a8 100644
--- a/src/math/remainder.c
+++ b/src/math/remainder.c
@@ -20,8 +20,6 @@
 
 #include "libm.h"
 
-static const double zero = 0.0;
-
 double remainder(double x, double p)
 {
 	int32_t hx,hp;
@@ -35,17 +33,15 @@ double remainder(double x, double p)
 	hx &= 0x7fffffff;
 
 	/* purge off exception values */
-	if ((hp|lp) == 0)  /* p = 0 */
-		return (x*p)/(x*p);
-	if (hx >= 0x7ff00000 ||                              /* x not finite */
+	if ((hp|lp) == 0 ||                                  /* p = 0 */
+	    hx >= 0x7ff00000 ||                              /* x not finite */
 	    (hp >= 0x7ff00000 && (hp-0x7ff00000 | lp) != 0)) /* p is NaN */
-		// FIXME: why long double?
-		return ((long double)x*p)/((long double)x*p);
+		return (x*p)/(x*p);
 
 	if (hp <= 0x7fdfffff)
 		x = fmod(x, p+p);  /* now x < 2p */
 	if (((hx-hp)|(lx-lp)) == 0)
-		return zero*x;
+		return 0.0*x;
 	x = fabs(x);
 	p = fabs(p);
 	if (hp < 0x00200000) {
diff --git a/src/math/remainderf.c b/src/math/remainderf.c
index 61c3c660..77671039 100644
--- a/src/math/remainderf.c
+++ b/src/math/remainderf.c
@@ -15,8 +15,6 @@
 
 #include "libm.h"
 
-static const float zero = 0.0;
-
 float remainderf(float x, float p)
 {
 	int32_t hx,hp;
@@ -30,16 +28,13 @@ float remainderf(float x, float p)
 	hx &= 0x7fffffff;
 
 	/* purge off exception values */
-	if (hp == 0)  /* p = 0 */
+	if (hp == 0 || hx >= 0x7f800000 || hp > 0x7f800000)  /* p = 0, x not finite, p is NaN */
 		return (x*p)/(x*p);
-	if (hx >= 0x7f800000 || hp > 0x7f800000)  /* x not finite, p is NaN */
-		// FIXME: why long double?
-		return ((long double)x*p)/((long double)x*p);
 
 	if (hp <= 0x7effffff)
 		x = fmodf(x, p + p);  /* now x < 2p */
 	if (hx - hp == 0)
-		return zero*x;
+		return 0.0f*x;
 	x = fabsf(x);
 	p = fabsf(p);
 	if (hp < 0x01000000) {
-- 
cgit v1.2.1


From 0cbb65479147ecdaa664e88cc2a5a925f3de502f Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Mon, 19 Mar 2012 23:41:19 +0100
Subject: code cleanup of named constants

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).
---
 src/math/__cos.c            |  5 ++--
 src/math/__cosdf.c          |  3 +--
 src/math/__cosl.c           |  6 ++---
 src/math/__rem_pio2.c       |  3 +--
 src/math/__rem_pio2_large.c | 12 ++++------
 src/math/__rem_pio2l.c      |  3 +--
 src/math/__sin.c            |  3 +--
 src/math/__sinl.c           |  4 +---
 src/math/__tan.c            | 10 +++-----
 src/math/acos.c             | 13 +++++------
 src/math/acosf.c            | 13 +++++------
 src/math/acosh.c            |  5 ++--
 src/math/acoshf.c           |  7 +++---
 src/math/acoshl.c           |  5 ++--
 src/math/acosl.c            |  5 ++--
 src/math/asin.c             |  9 ++++---
 src/math/asinf.c            |  9 ++++---
 src/math/asinh.c            |  7 +++---
 src/math/asinhf.c           |  7 +++---
 src/math/asinhl.c           |  7 +++---
 src/math/asinl.c            |  8 +++----
 src/math/atan.c             | 12 ++++------
 src/math/atan2.c            |  5 ++--
 src/math/atan2f.c           |  5 ++--
 src/math/atan2l.c           |  5 ++--
 src/math/atanf.c            | 12 ++++------
 src/math/atanh.c            | 11 ++++-----
 src/math/atanhf.c           | 11 ++++-----
 src/math/atanhl.c           | 10 ++++----
 src/math/atanl.c            | 12 ++++------
 src/math/cosh.c             | 12 +++++-----
 src/math/coshf.c            | 14 +++++------
 src/math/coshl.c            | 16 ++++++-------
 src/math/erf.c              | 57 +++++++++++++++++++++------------------------
 src/math/erff.c             | 57 +++++++++++++++++++++------------------------
 src/math/erfl.c             | 54 ++++++++++++++++++++----------------------
 src/math/exp.c              |  9 ++++---
 src/math/expf.c             |  9 ++++---
 src/math/expm1.c            | 15 ++++++------
 src/math/expm1f.c           | 15 ++++++------
 src/math/j0.c               | 39 ++++++++++++++-----------------
 src/math/j0f.c              | 39 ++++++++++++++-----------------
 src/math/j1.c               | 33 ++++++++++++--------------
 src/math/j1f.c              | 33 ++++++++++++--------------
 src/math/jn.c               | 33 +++++++++++---------------
 src/math/jnf.c              | 32 +++++++++++--------------
 src/math/lgamma_r.c         | 50 ++++++++++++++++++---------------------
 src/math/lgammaf_r.c        | 40 ++++++++++++++-----------------
 src/math/lgammal.c          | 48 +++++++++++++++++---------------------
 src/math/log.c              | 10 ++++----
 src/math/log10.c            |  8 +++----
 src/math/log10f.c           |  8 +++----
 src/math/log1p.c            | 10 ++++----
 src/math/log1pf.c           | 10 ++++----
 src/math/log1pl.c           | 20 ++++++++--------
 src/math/log2.c             |  8 +++----
 src/math/log2f.c            |  8 +++----
 src/math/logf.c             | 10 ++++----
 src/math/modf.c             |  4 +---
 src/math/pow.c              | 37 ++++++++++++++---------------
 src/math/powf.c             | 31 +++++++++++-------------
 src/math/remquol.c          |  2 +-
 src/math/rintl.c            |  2 +-
 src/math/scalbnl.c          |  2 +-
 src/math/sinh.c             |  8 +++----
 src/math/sinhf.c            |  8 +++----
 src/math/sinhl.c            |  8 +++----
 src/math/sqrt.c             | 10 ++++----
 src/math/sqrtf.c            | 10 ++++----
 src/math/tanh.c             | 18 +++++++-------
 src/math/tanhf.c            | 20 +++++++++-------
 src/math/tanhl.c            | 18 +++++++-------
 src/math/tgammal.c          | 34 +++++++++++++--------------
 73 files changed, 513 insertions(+), 623 deletions(-)

(limited to 'src/math')

diff --git a/src/math/__cos.c b/src/math/__cos.c
index ba439857..8699c1d5 100644
--- a/src/math/__cos.c
+++ b/src/math/__cos.c
@@ -51,7 +51,6 @@
 #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 */
@@ -67,6 +66,6 @@ double __cos(double x, double y)
 	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));
 }
diff --git a/src/math/__cosdf.c b/src/math/__cosdf.c
index a3b399e6..a65f7f21 100644
--- a/src/math/__cosdf.c
+++ b/src/math/__cosdf.c
@@ -18,7 +18,6 @@
 
 /* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
 static const double
-one =  1.0,
 C0  = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
 C1  =  0x155553e1053a42.0p-57, /*  0.0416666233237390631894 */
 C2  = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
@@ -32,5 +31,5 @@ float __cosdf(double x)
 	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;
 }
diff --git a/src/math/__cosl.c b/src/math/__cosl.c
index 80036ddb..9d325768 100644
--- a/src/math/__cosl.c
+++ b/src/math/__cosl.c
@@ -41,8 +41,6 @@
  * 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 */
 
@@ -61,7 +59,7 @@ long double __cosl(long double x, long double y)
 	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));
 }
 #endif
diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c
index a7d779e0..0ef57fb5 100644
--- a/src/math/__rem_pio2.c
+++ b/src/math/__rem_pio2.c
@@ -29,7 +29,6 @@
  * 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 */
@@ -163,7 +162,7 @@ medium:
 	}
 	tx[2] = z;
 	nx = 3;
-	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];
diff --git a/src/math/__rem_pio2_large.c b/src/math/__rem_pio2_large.c
index 35835f83..27b619cc 100644
--- a/src/math/__rem_pio2_large.c
+++ b/src/math/__rem_pio2_large.c
@@ -271,8 +271,6 @@ static const double PIo2[] = {
 };
 
 static const double
-zero   = 0.0,
-one    = 1.0,
 two24  = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
 
@@ -293,7 +291,7 @@ int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
 	/* 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++) {
@@ -346,14 +344,14 @@ recompute:
 			}
 		}
 		if (ih == 2) {
-			z = one - z;
+			z = 1.0 - z;
 			if (carry != 0)
-				z -= scalbn(one,q0);
+				z -= scalbn(1.0,q0);
 		}
 	}
 
 	/* check if recomputation is needed */
-	if (z == zero) {
+	if (z == 0.0) {
 		j = 0;
 		for (i=jz-1; i>=jk; i--) j |= iq[i];
 		if (j == 0) {  /* need recomputation */
@@ -391,7 +389,7 @@ recompute:
 	}
 
 	/* convert integer "bit" chunk to floating-point value */
-	fw = scalbn(one,q0);
+	fw = scalbn(1.0,q0);
 	for (i=jz; i>=0; i--) {
 		q[i] = fw*(double)iq[i];
 		fw *= twon24;
diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c
index 10af404c..f0cb99ac 100644
--- a/src/math/__rem_pio2l.c
+++ b/src/math/__rem_pio2l.c
@@ -32,7 +32,6 @@
  * 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 */
@@ -114,7 +113,7 @@ int __rem_pio2l(long double x, long double *y)
 	}
 	tx[2] = z;
 	nx = 3;
-	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];
diff --git a/src/math/__sin.c b/src/math/__sin.c
index 80f3273c..9aead04b 100644
--- a/src/math/__sin.c
+++ b/src/math/__sin.c
@@ -42,7 +42,6 @@
 #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 */
@@ -61,5 +60,5 @@ double __sin(double x, double y, int iy)
 	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);
 }
diff --git a/src/math/__sinl.c b/src/math/__sinl.c
index 67c4bdc5..068adffb 100644
--- a/src/math/__sinl.c
+++ b/src/math/__sinl.c
@@ -24,8 +24,6 @@
  * 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 */
 
@@ -47,6 +45,6 @@ long double __sinl(long double x, long double y, int iy)
 	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);
 }
 #endif
diff --git a/src/math/__tan.c b/src/math/__tan.c
index f1be2ec8..01e3fe48 100644
--- a/src/math/__tan.c
+++ b/src/math/__tan.c
@@ -59,13 +59,9 @@ static const double T[] = {
              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)
 {
diff --git a/src/math/acos.c b/src/math/acos.c
index 456a2219..54d266ee 100644
--- a/src/math/acos.c
+++ b/src/math/acos.c
@@ -36,8 +36,7 @@
 #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 */
@@ -75,25 +74,25 @@ double acos(double x)
 			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 */
-		z = (one+x)*0.5;
+		z = (1.0+x)*0.5;
 		p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
-		q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+		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 */
-		z = (one-x)*0.5;
+		z = (1.0-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)))));
-		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);
diff --git a/src/math/acosf.c b/src/math/acosf.c
index 945347cb..d875f3ef 100644
--- a/src/math/acosf.c
+++ b/src/math/acosf.c
@@ -16,8 +16,7 @@
 #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 */
@@ -46,13 +45,13 @@ float acosf(float x)
 			return pio2_hi + pio2_lo;
 		z = x*x;
 		p = z*(pS0+z*(pS1+z*pS2));
-		q = one+z*qS1;
+		q = 1.0f+z*qS1;
 		r = p/q;
 		return pio2_hi - (x - (pio2_lo-x*r));
 	} else if (hx < 0) {     /* x < -0.5 */
-		z = (one+x)*0.5f;
+		z = (1.0f+x)*0.5f;
 		p = z*(pS0+z*(pS1+z*pS2));
-		q = one+z*qS1;
+		q = 1.0f+z*qS1;
 		s = sqrtf(z);
 		r = p/q;
 		w = r*s-pio2_lo;
@@ -60,14 +59,14 @@ float acosf(float x)
 	} else {                 /* x > 0.5 */
 		int32_t idf;
 
-		z = (one-x)*0.5f;
+		z = (1.0f-x)*0.5f;
 		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));
-		q = one+z*qS1;
+		q = 1.0f+z*qS1;
 		r = p/q;
 		w = r*s+c;
 		return 2.0f*(df+w);
diff --git a/src/math/acosh.c b/src/math/acosh.c
index a7c87e3c..15f51c6e 100644
--- a/src/math/acosh.c
+++ b/src/math/acosh.c
@@ -27,7 +27,6 @@
 #include "libm.h"
 
 static const double
-one = 1.0,
 ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
 
 double acosh(double x)
@@ -47,9 +46,9 @@ double acosh(double 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)));
 	} else {                /* 1 < x < 2 */
-		t = x-one;
+		t = x-1.0;
 		return log1p(t + sqrt(2.0*t+t*t));
 	}
 }
diff --git a/src/math/acoshf.c b/src/math/acoshf.c
index 57ce5cb8..0f7aae2a 100644
--- a/src/math/acoshf.c
+++ b/src/math/acoshf.c
@@ -16,7 +16,6 @@
 #include "libm.h"
 
 static const float
-one = 1.0,
 ln2 = 6.9314718246e-01; /* 0x3f317218 */
 
 float acoshf(float x)
@@ -32,12 +31,12 @@ float acoshf(float x)
 			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;
-		return logf(2.0f*x - one/(x+sqrtf(t-one)));
+		return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f)));
 	} else {                /* 1 < x < 2 */
-		t = x-one;
+		t = x-1.0f;
 		return log1pf(t + sqrtf(2.0f*t+t*t));
 	}
 }
diff --git a/src/math/acoshl.c b/src/math/acoshl.c
index d8310a73..a4024516 100644
--- a/src/math/acoshl.c
+++ b/src/math/acoshl.c
@@ -32,7 +32,6 @@ long double acoshl(long double x)
 }
 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
 static const long double
-one = 1.0,
 ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
 
 long double acoshl(long double x)
@@ -51,10 +50,10 @@ long double acoshl(long double 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)));
 	}
 	/* 1 < x <= 2 */
-	t = x - one;
+	t = x - 1.0;
 	return log1pl(t + sqrtl(2.0*t + t*t));
 }
 #endif
diff --git a/src/math/acosl.c b/src/math/acosl.c
index 170520fe..cc565336 100644
--- a/src/math/acosl.c
+++ b/src/math/acosl.c
@@ -25,7 +25,6 @@ long double acosl(long double x)
 #include "__invtrigl.h"
 
 static const long double
-one = 1.00000000000000000000e+00,
 pi = 3.14159265358979323846264338327950280e+00L;
 
 long double acosl(long double x)
@@ -55,7 +54,7 @@ long double acosl(long double x)
 		r = p / q;
 		return pio2_hi - (x - (pio2_lo - x * r));
 	} else if (expsign < 0) {  /* x < -0.5 */
-		z = (one + x) * 0.5;
+		z = (1.0 + x) * 0.5;
 		p = P(z);
 		q = Q(z);
 		s = sqrtl(z);
@@ -63,7 +62,7 @@ long double acosl(long double x)
 		w = r * s - pio2_lo;
 		return pi - 2.0 * (s + w);
 	} else {                   /* x > 0.5 */
-		z = (one - x) * 0.5;
+		z = (1.0 - x) * 0.5;
 		s = sqrtl(z);
 		u.e = s;
 		u.bits.manl = 0;
diff --git a/src/math/asin.c b/src/math/asin.c
index 04bd0c14..e3d8b250 100644
--- a/src/math/asin.c
+++ b/src/math/asin.c
@@ -42,7 +42,6 @@
 #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 */
@@ -76,20 +75,20 @@ double asin(double x)
 		return (x-x)/(x-x);  /* asin(|x|>1) is NaN */
 	} else if (ix < 0x3fe00000) {  /* |x|<0.5 */
 		if (ix < 0x3e500000) {  /* if |x| < 2**-26 */
-			if (huge+x > one)
+			if (huge+x > 1.0)
 				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 = one - fabs(x);
+	w = 1.0 - fabs(x);
 	t = w*0.5;
 	p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
-	q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+	q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
 	s = sqrt(t);
 	if (ix >= 0x3FEF3333) {  /* if |x| > 0.975 */
 		w = p/q;
diff --git a/src/math/asinf.c b/src/math/asinf.c
index 7865a45b..c1c42c7b 100644
--- a/src/math/asinf.c
+++ b/src/math/asinf.c
@@ -16,7 +16,6 @@
 #include "libm.h"
 
 static const float
-one =  1.0000000000e+00, /* 0x3F800000 */
 huge = 1.000e+30,
 /* coefficients for R(x^2) */
 pS0 =  1.6666586697e-01,
@@ -41,20 +40,20 @@ float asinf(float x)
 		return (x-x)/(x-x);  /* asin(|x|>1) is NaN */
 	} else if (ix < 0x3f000000) {  /* |x|<0.5 */
 		if (ix < 0x39800000) {  /* |x| < 2**-12 */
-			if (huge+x > one)
+			if (huge+x > 1.0f)
 				return x; /* return x with inexact if x!=0 */
 		}
 		t = x*x;
 		p = t*(pS0+t*(pS1+t*pS2));
-		q = one+t*qS1;
+		q = 1.0f+t*qS1;
 		w = p/q;
 		return x + x*w;
 	}
 	/* 1 > |x| >= 0.5 */
-	w = one - fabsf(x);
+	w = 1.0f - fabsf(x);
 	t = w*0.5f;
 	p = t*(pS0+t*(pS1+t*pS2));
-	q = one+t*qS1;
+	q = 1.0f+t*qS1;
 	s = sqrt(t);
 	w = p/q;
 	t = pio2-2.0*(s+s*w);
diff --git a/src/math/asinh.c b/src/math/asinh.c
index 92aa9446..11bbd71a 100644
--- a/src/math/asinh.c
+++ b/src/math/asinh.c
@@ -23,7 +23,6 @@
 #include "libm.h"
 
 static const double
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
 huge= 1.00000000000000000000e+300;
 
@@ -38,17 +37,17 @@ double asinh(double x)
 		return x+x;
 	if (ix < 0x3e300000) {  /* |x| < 2**-28 */
 		/* return x inexact except 0 */
-		if (huge+x > one)
+		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 + 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;
-		w =log1p(fabs(x) + t/(one+sqrt(one+t)));
+		w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
 	}
 	if (hx > 0)
 		return w;
diff --git a/src/math/asinhf.c b/src/math/asinhf.c
index cb6e5b89..efe3af94 100644
--- a/src/math/asinhf.c
+++ b/src/math/asinhf.c
@@ -16,7 +16,6 @@
 #include "libm.h"
 
 static const float
-one = 1.0000000000e+00, /* 0x3F800000 */
 ln2 = 6.9314718246e-01, /* 0x3f317218 */
 huge= 1.0000000000e+30;
 
@@ -31,17 +30,17 @@ float asinhf(float x)
 		return x+x;
 	if (ix < 0x31800000) {  /* |x| < 2**-28 */
 		/* return x inexact except 0 */
-		if (huge+x > one)
+		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 + 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;
-		w =log1pf(fabsf(x) + t/(one+sqrtf(one+t)));
+		w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
 	}
 	if (hx > 0)
 		return w;
diff --git a/src/math/asinhl.c b/src/math/asinhl.c
index b2edf904..dc5dd71f 100644
--- a/src/math/asinhl.c
+++ b/src/math/asinhl.c
@@ -29,7 +29,6 @@ long double asinhl(long double x)
 }
 #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;
 
@@ -44,17 +43,17 @@ long double asinhl(long double x)
 		return x + x;   /* x is inf or NaN */
 	if (ix < 0x3fde) {      /* |x| < 2**-34 */
 		/* return x, raise inexact if x != 0 */
-		if (huge+x > one)
+		if (huge+x > 1.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);
-		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;
-		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;
diff --git a/src/math/asinl.c b/src/math/asinl.c
index 370997bc..ddd807e2 100644
--- a/src/math/asinl.c
+++ b/src/math/asinl.c
@@ -23,9 +23,7 @@ long double asinl(long double x)
 }
 #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)
 {
@@ -45,7 +43,7 @@ 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 */
-			if (huge+x > one)
+			if (huge+x > 1.0)
 				return x;
 		}
 		t = x*x;
@@ -55,7 +53,7 @@ long double asinl(long double x)
 		return x + x*w;
 	}
 	/* 1 > |x| >= 0.5 */
-	w = one - fabsl(x);
+	w = 1.0 - fabsl(x);
 	t = w*0.5;
 	p = P(t);
 	q = Q(t);
diff --git a/src/math/atan.c b/src/math/atan.c
index d31782c2..f34551c7 100644
--- a/src/math/atan.c
+++ b/src/math/atan.c
@@ -60,9 +60,7 @@ static const double aT[] = {
   1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
 };
 
-static const double
-one = 1.0,
-huge = 1.0e300;
+static const double huge = 1.0e300;
 
 double atan(double x)
 {
@@ -86,7 +84,7 @@ double atan(double x)
 	if (ix < 0x3fdc0000) {    /* |x| < 0.4375 */
 		if (ix < 0x3e400000) {  /* |x| < 2^-27 */
 			/* raise inexact */
-			if (huge+x > one)
+			if (huge+x > 1.0)
 				return x;
 		}
 		id = -1;
@@ -95,15 +93,15 @@ double atan(double x)
 		if (ix < 0x3ff30000) {  /* |x| < 1.1875 */
 			if (ix < 0x3fe60000) {  /*  7/16 <= |x| < 11/16 */
 				id = 0;
-				x = (2.0*x-one)/(2.0+x);
+				x = (2.0*x-1.0)/(2.0+x);
 			} else {                /* 11/16 <= |x| < 19/16 */
 				id = 1;
-				x = (x-one)/(x+one);
+				x = (x-1.0)/(x+1.0);
 			}
 		} 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;
diff --git a/src/math/atan2.c b/src/math/atan2.c
index d928f0ed..143c3834 100644
--- a/src/math/atan2.c
+++ b/src/math/atan2.c
@@ -42,7 +42,6 @@
 static const volatile double
 tiny  = 1.0e-300;
 static const double
-zero  = 0.0,
 pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
 pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
 pi     = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
@@ -89,8 +88,8 @@ double atan2(double y, double x)
 			}
 		} else {
 			switch(m) {
-			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) */
 			}
diff --git a/src/math/atan2f.c b/src/math/atan2f.c
index 19b134dc..96839cba 100644
--- a/src/math/atan2f.c
+++ b/src/math/atan2f.c
@@ -18,7 +18,6 @@
 static const volatile float
 tiny = 1.0e-30;
 static const float
-zero = 0.0,
 pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
 pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
 pi     = 3.1415927410e+00; /* 0x40490fdb */
@@ -63,8 +62,8 @@ float atan2f(float y, float x)
 			}
 		} else {
 			switch (m) {
-			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) */
 			}
diff --git a/src/math/atan2l.c b/src/math/atan2l.c
index 48abc058..0fc901c8 100644
--- a/src/math/atan2l.c
+++ b/src/math/atan2l.c
@@ -27,7 +27,6 @@ long double atan2l(long double y, long double x)
 static const volatile long double
 tiny = 1.0e-300;
 static const long double
-zero = 0.0,
 pi = 3.14159265358979323846264338327950280e+00L;
 
 long double atan2l(long double y, long double x)
@@ -75,8 +74,8 @@ long double atan2l(long double y, long double x)
 			}
 		} else {
 			switch(m) {
-			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) */
 			}
diff --git a/src/math/atanf.c b/src/math/atanf.c
index 73cebb5e..4e31b902 100644
--- a/src/math/atanf.c
+++ b/src/math/atanf.c
@@ -38,9 +38,7 @@ static const float aT[] = {
   6.1687607318e-02,
 };
 
-static const float
-one = 1.0,
-huge = 1.0e30;
+static const float huge = 1.0e30;
 
 float atanf(float x)
 {
@@ -60,7 +58,7 @@ float atanf(float x)
 	if (ix < 0x3ee00000) {   /* |x| < 0.4375 */
 		if (ix < 0x39800000) {  /* |x| < 2**-12 */
 			/* raise inexact */
-			if(huge+x>one)
+			if(huge+x>1.0f)
 				return x;
 		}
 		id = -1;
@@ -69,15 +67,15 @@ float atanf(float x)
 		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;
-				x = (x - one)/(x + one);
+				x = (x - 1.0f)/(x + 1.0f);
 			}
 		} 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;
diff --git a/src/math/atanh.c b/src/math/atanh.c
index 29290463..dbe241d1 100644
--- a/src/math/atanh.c
+++ b/src/math/atanh.c
@@ -30,8 +30,7 @@
 
 #include "libm.h"
 
-static const double one = 1.0, huge = 1e300;
-static const double zero = 0.0;
+static const double huge = 1e300;
 
 double atanh(double x)
 {
@@ -44,15 +43,15 @@ double atanh(double x)
 	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;
-		t = 0.5*log1p(t + t*x/(one-x));
+		t = 0.5*log1p(t + t*x/(1.0-x));
 	} else
-		t = 0.5*log1p((x+x)/(one-x));
+		t = 0.5*log1p((x+x)/(1.0-x));
 	if (hx >= 0)
 		return t;
 	return -t;
diff --git a/src/math/atanhf.c b/src/math/atanhf.c
index 9c65dc7f..2be780bb 100644
--- a/src/math/atanhf.c
+++ b/src/math/atanhf.c
@@ -15,8 +15,7 @@
 
 #include "libm.h"
 
-static const float one = 1.0, huge = 1e30;
-static const float zero = 0.0;
+static const float huge = 1e30;
 
 float atanhf(float x)
 {
@@ -28,15 +27,15 @@ float atanhf(float x)
 	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;
-		t = 0.5f*log1pf(t + t*x/(one-x));
+		t = 0.5f*log1pf(t + t*x/(1.0f-x));
 	} else
-		t = 0.5f*log1pf((x+x)/(one-x));
+		t = 0.5f*log1pf((x+x)/(1.0f-x));
 	if (hx >= 0)
 		return t;
 	return -t;
diff --git a/src/math/atanhl.c b/src/math/atanhl.c
index af0f856d..931bae32 100644
--- a/src/math/atanhl.c
+++ b/src/math/atanhl.c
@@ -34,7 +34,7 @@ long double atanhl(long double x)
 	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)
 {
@@ -48,15 +48,15 @@ long double atanhl(long double x)
 		/* |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;
-		t = 0.5*log1pl(t + t*x/(one-x));
+		t = 0.5*log1pl(t + t*x/(1.0 - x));
 	} else
-		t = 0.5*log1pl((x + x)/(one - x));
+		t = 0.5*log1pl((x + x)/(1.0 - x));
 	if (se <= 0x7fff)
 		return t;
 	return -t;
diff --git a/src/math/atanl.c b/src/math/atanl.c
index 4e99955e..36072c17 100644
--- a/src/math/atanl.c
+++ b/src/math/atanl.c
@@ -23,9 +23,7 @@ long double atanl(long double x)
 }
 #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)
 {
@@ -53,7 +51,7 @@ 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 (huge+x > one)
+			if (huge+x > 1.0)
 				return x;
 		}
 		id = -1;
@@ -62,15 +60,15 @@ long double atanl(long double x)
 		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;
-				x = (x-one)/(x+one);
+				x = (x-1.0)/(x+1.0);
 			}
 		} 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;
diff --git a/src/math/cosh.c b/src/math/cosh.c
index 5f38b276..a1f7dbc7 100644
--- a/src/math/cosh.c
+++ b/src/math/cosh.c
@@ -32,7 +32,7 @@
 
 #include "libm.h"
 
-static const double one = 1.0, half = 0.5, huge = 1.0e300;
+static const double huge = 1.0e300;
 
 double cosh(double x)
 {
@@ -49,21 +49,21 @@ double cosh(double x)
 	/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
 	if (ix < 0x3fd62e43) {
 		t = expm1(fabs(x));
-		w = one+t;
+		w = 1.0+t;
 		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));
-		return half*t + half/t;
+		return 0.5*t + 0.5/t;
 	}
 
-	/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+	/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
 	if (ix < 0x40862E42)
-		return half*exp(fabs(x));
+		return 0.5*exp(fabs(x));
 
 	/* |x| in [log(maxdouble), overflowthresold] */
 	if (ix <= 0x408633CE)
diff --git a/src/math/coshf.c b/src/math/coshf.c
index 9e87afcd..97318f10 100644
--- a/src/math/coshf.c
+++ b/src/math/coshf.c
@@ -15,7 +15,7 @@
 
 #include "libm.h"
 
-static const float one = 1.0, half = 0.5, huge = 1.0e30;
+static const float huge = 1.0e30;
 
 float coshf(float x)
 {
@@ -32,21 +32,21 @@ float coshf(float x)
 	/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
 	if (ix < 0x3eb17218) {
 		t = expm1f(fabsf(x));
-		w = one+t;
+		w = 1.0f+t;
 		if (ix<0x39800000)
-			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));
-		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|) */
 	if (ix < 0x42b17217)
-		return half*expf(fabsf(x));
+		return 0.5f*expf(fabsf(x));
 
 	/* |x| in [log(maxfloat), overflowthresold] */
 	if (ix <= 0x42b2d4fc)
diff --git a/src/math/coshl.c b/src/math/coshl.c
index bcc9128a..36b52c02 100644
--- a/src/math/coshl.c
+++ b/src/math/coshl.c
@@ -38,7 +38,7 @@ long double coshl(long double x)
 	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)
 {
@@ -56,27 +56,27 @@ 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));
-		w = one + t;
+		w = 1.0 + t;
 		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));
-		return half*t + half/t;
+		return 0.5*t + 0.5/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))
-		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)))
 	{
-		w = expl(half*fabsl(x));
-		t = half*w;
+		w = expl(0.5*fabsl(x));
+		t = 0.5*w;
 		return t*w;
 	}
 
diff --git a/src/math/erf.c b/src/math/erf.c
index 18ee01cf..dd8c3a77 100644
--- a/src/math/erf.c
+++ b/src/math/erf.c
@@ -107,9 +107,6 @@
 
 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 */
 /*
@@ -190,7 +187,7 @@ double erf(double x)
 	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 */
@@ -201,42 +198,42 @@ double erf(double x)
 		}
 		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 */
-		s = fabs(x)-one;
+		s = fabs(x)-1.0;
 		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)
-			return one-tiny;
-		return tiny-one;
+			return 1.0 - tiny;
+		return tiny - 1.0;
 	}
 	x = fabs(x);
-	s = one/(x*x);
+	s = 1.0/(x*x);
 	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)))));
-		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)
-		return one-r/x;
-	return r/x-one;
+		return 1.0 - r/x;
+	return r/x - 1.0;
 }
 
 double erfc(double x)
@@ -248,49 +245,49 @@ double erfc(double x)
 	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 */
-			return one - x;
+			return 1.0 - x;
 		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 */
-			return one - (x+x*y);
+			return 1.0 - (x+x*y);
 		} else {
 			r = x*y;
-			r += x-half;
-			return half - r ;
+			r += x - 0.5;
+			return 0.5 - r ;
 		}
 	}
 	if (ix < 0x3ff40000) {  /* 0.84375 <= |x| < 1.25 */
-		s = fabs(x)-one;
+		s = fabs(x)-1.0;
 		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) {
-			z = one-erx;
+			z = 1.0-erx;
 			return z - P/Q;
 		} else {
 			z = erx+P/Q;
-			return one+z;
+			return 1.0 + z;
 		}
 	}
 	if (ix < 0x403c0000) {  /* |x| < 28 */
 		x = fabs(x);
-		s = one/(x*x);
+		s = 1.0/(x*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))))));
-			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 */
-				return two-tiny;
+				return 2.0 - tiny;
 			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;
@@ -298,9 +295,9 @@ double erfc(double x)
 		r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 		if (hx > 0)
 			return r/x;
-		return two-r/x;
+		return 2.0 - r/x;
 	}
 	if (hx > 0)
 		return tiny*tiny;
-	return two-tiny;
+	return 2.0 - tiny;
 }
diff --git a/src/math/erff.c b/src/math/erff.c
index eef4851e..a0a304e5 100644
--- a/src/math/erff.c
+++ b/src/math/erff.c
@@ -17,9 +17,6 @@
 
 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 */
 /*
@@ -100,7 +97,7 @@ float erff(float x)
 	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 */
@@ -111,42 +108,42 @@ float erff(float x)
 		}
 		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 */
-		s = fabsf(x)-one;
+		s = fabsf(x)-1.0f;
 		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)
-			return one - tiny;
-		return tiny - one;
+			return 1.0f - tiny;
+		return tiny - 1.0f;
 	}
 	x = fabsf(x);
-	s = one/(x*x);
+	s = 1.0f/(x*x);
 	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)))));
-		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)
-		return one - r/x;
-	return  r/x - one;
+		return 1.0f - r/x;
+	return  r/x - 1.0f;
 }
 
 float erfcf(float x)
@@ -158,50 +155,50 @@ float erfcf(float x)
 	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 */
-			return one - x;
+			return 1.0f - x;
 		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 */
-			return one - (x+x*y);
+			return 1.0f - (x+x*y);
 		} else {
 			r = x*y;
-			r += (x-half);
-			return half - r ;
+			r += (x-0.5f);
+			return 0.5f - r ;
 		}
 	}
 	if (ix < 0x3fa00000) {  /* 0.84375 <= |x| < 1.25 */
-		s = fabsf(x)-one;
+		s = fabsf(x)-1.0f;
 		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) {
-			z = one - erx;
+			z = 1.0f - erx;
 			return z - P/Q;
 		} else {
 			z = erx + P/Q;
-			return one + z;
+			return 1.0f + z;
 		}
 	}
 	if (ix < 0x41e00000) {  /* |x| < 28 */
 		x = fabsf(x);
-		s = one/(x*x);
+		s = 1.0f/(x*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))))));
-			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 */
-				return two-tiny;
+				return 2.0f-tiny;
 			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);
@@ -209,9 +206,9 @@ float erfcf(float x)
 		r = expf(-z*z - 0.5625f) * expf((z-x)*(z+x) + R/S);
 		if (hx > 0)
 			return r/x;
-		return two - r/x;
+		return 2.0f - r/x;
 	}
 	if (hx > 0)
 		return tiny*tiny;
-	return two - tiny;
+	return 2.0f - tiny;
 }
diff --git a/src/math/erfl.c b/src/math/erfl.c
index a80c2ce1..f7449859 100644
--- a/src/math/erfl.c
+++ b/src/math/erfl.c
@@ -108,9 +108,6 @@ long double erfl(long double x)
 #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,
 
@@ -248,12 +245,12 @@ long double erfl(long double x)
 	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;
-		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);
@@ -272,7 +269,7 @@ long double erfl(long double x)
 		return x + x * y;
 	}
 	if (ix < 0x3fffa000) {  /* 0.84375 <= |x| < 1.25 */
-		s = fabsl (x) - one;
+		s = fabsl(x) - 1.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] +
@@ -283,11 +280,11 @@ long double erfl(long double x)
 	}
 	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;
 	}
 	x = fabsl (x);
-	s = one / (x * 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])))))));
@@ -305,8 +302,8 @@ long double erfl(long double x)
 	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)
@@ -315,16 +312,16 @@ long double erfcl(long double x)
 	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 */
-		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 */
-			return one - x;
+			return 1.0 - x;
 		z = x * x;
 		r = pp[0] + z * (pp[1] +
 		     z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
@@ -332,27 +329,27 @@ long double erfcl(long double x)
 		     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 * y;
-		r += x - half;
-		return half - r;
+		r += x - 0.5L;
+		return 0.5L - r;
 	}
 	if (ix < 0x3fffa000) {  /* 0.84375 <= |x| < 1.25 */
-		s = fabsl (x) - one;
+		s = fabsl(x) - 1.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) {
-			z = one - erx;
+			z = 1.0 - erx;
 			return z - P / Q;
 		}
 		z = erx + P / Q;
-		return one + z;
+		return 1.0 + z;
 	}
 	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])))))));
@@ -365,26 +362,25 @@ long double erfcl(long double x)
 			     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;
-		GET_LDOUBLE_WORDS (hx, i0, i1, z);
+		GET_LDOUBLE_WORDS(hx, i0, i1, z);
 		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;
-		return two - r / x;
+		return 2.0 - r / x;
 	}
 
 	if ((se & 0x8000) == 0)
 		return tiny * tiny;
-	return two - tiny;
+	return 2.0 - tiny;
 }
 #endif
diff --git a/src/math/exp.c b/src/math/exp.c
index a538b8cd..29bf9609 100644
--- a/src/math/exp.c
+++ b/src/math/exp.c
@@ -74,7 +74,6 @@
 #include "libm.h"
 
 static const double
-one     = 1.0,
 halF[2] = {0.5,-0.5,},
 huge    = 1.0e+300,
 o_threshold =  7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
@@ -134,8 +133,8 @@ double exp(double x)
 		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;
 	} else
 		k = 0;
 
@@ -147,8 +146,8 @@ double exp(double x)
 		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)
diff --git a/src/math/expf.c b/src/math/expf.c
index f706ac5d..3c3e2ab9 100644
--- a/src/math/expf.c
+++ b/src/math/expf.c
@@ -16,7 +16,6 @@
 #include "libm.h"
 
 static const float
-one     = 1.0,
 halF[2] = {0.5,-0.5,},
 huge    = 1.0e+30,
 o_threshold =  8.8721679688e+01,  /* 0x42b17180 */
@@ -72,8 +71,8 @@ float expf(float x)
 		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;
 	} else
 		k = 0;
 
@@ -85,8 +84,8 @@ float expf(float x)
 		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)
diff --git a/src/math/expm1.c b/src/math/expm1.c
index ffa82264..f8f32c46 100644
--- a/src/math/expm1.c
+++ b/src/math/expm1.c
@@ -107,7 +107,6 @@
 #include "libm.h"
 
 static const double
-one         = 1.0,
 huge        = 1.0e+300,
 tiny        = 1.0e-300,
 o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
@@ -148,7 +147,7 @@ double expm1(double x)
 		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 */
 		}
 	}
 
@@ -182,7 +181,7 @@ double expm1(double x)
 	/* 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 */
@@ -195,17 +194,17 @@ double expm1(double x)
 	if (k == 1) {
 		if (x < -0.25)
 			return -2.0*(e-(x+0.5));
-		return one+2.0*(x-e);
+		return 1.0+2.0*(x-e);
 	}
 	if (k <= -2 || k > 56) {  /* suffice to return exp(x)-1 */
-		y = one - (e-x);
+		y = 1.0 - (e-x);
 		if (k == 1024)
 			y = y*2.0*0x1p1023;
 		else
 			y = y*twopk;
-		return y - one;
+		return y - 1.0;
 	}
-	t = one;
+	t = 1.0;
 	if (k < 20) {
 		SET_HIGH_WORD(t, 0x3ff00000 - (0x200000>>k));  /* t=1-2^-k */
 		y = t-(e-x);
@@ -213,7 +212,7 @@ double expm1(double x)
 	} else {
 		SET_HIGH_WORD(t, ((0x3ff-k)<<20));  /* 2^-k */
 		y = x-(e+t);
-		y += one;
+		y += 1.0;
 		y = y*twopk;
 	}
 	return y;
diff --git a/src/math/expm1f.c b/src/math/expm1f.c
index a8b79e88..d9568466 100644
--- a/src/math/expm1f.c
+++ b/src/math/expm1f.c
@@ -16,7 +16,6 @@
 #include "libm.h"
 
 static const float
-one         = 1.0,
 huge        = 1.0e+30,
 tiny        = 1.0e-30,
 o_threshold = 8.8721679688e+01, /* 0x42b17180 */
@@ -54,7 +53,7 @@ float expm1f(float x)
 		if (xsb != 0) {  /* x < -27*ln2 */
 			/* raise inexact */
 			if (x+tiny < 0.0f)
-				return tiny-one;  /* return -1 */
+				return tiny-1.0f;  /* return -1 */
 		}
 	}
 
@@ -87,7 +86,7 @@ float expm1f(float x)
 	/* 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 */
@@ -100,17 +99,17 @@ float expm1f(float x)
 	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 */
-		y = one - (e - x);
+		y = 1.0f - (e - x);
 		if (k == 128)
 			y = y*2.0f*0x1p127f;
 		else
 			y = y*twopk;
-		return y - one;
+		return y - 1.0f;
 	}
-	t = one;
+	t = 1.0f;
 	if (k < 23) {
 		SET_FLOAT_WORD(t, 0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
 		y = t - (e - x);
@@ -118,7 +117,7 @@ float expm1f(float x)
 	} else {
 		SET_FLOAT_WORD(t, (0x7f-k)<<23);  /* 2^-k */
 		y = x - (e + t);
-		y += one;
+		y += 1.0f;
 		y = y*twopk;
 	}
 	return y;
diff --git a/src/math/j0.c b/src/math/j0.c
index b5490641..986968e8 100644
--- a/src/math/j0.c
+++ b/src/math/j0.c
@@ -60,7 +60,6 @@ static double pzero(double), qzero(double);
 
 static const double
 huge      = 1e300,
-one       = 1.0,
 invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
 tpi       = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 /* R0/S0 on [0, 2.00] */
@@ -73,8 +72,6 @@ S02 =  1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
 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;
@@ -83,7 +80,7 @@ double j0(double x)
 	GET_HIGH_WORD(hx, x);
 	ix = hx & 0x7fffffff;
 	if (ix >= 0x7ff00000)
-		return one/(x*x);
+		return 1.0/(x*x);
 	x = fabs(x);
 	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
 		s = sin(x);
@@ -92,7 +89,7 @@ double j0(double x)
 		cc = s+c;
 		if (ix < 0x7fe00000) {  /* make sure x+x does not overflow */
 			z = -cos(x+x);
-			if ((s*c) < zero)
+			if (s*c < 0.0)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -112,20 +109,20 @@ double j0(double x)
 	}
 	if (ix < 0x3f200000) {  /* |x| < 2**-13 */
 		/* raise inexact if x != 0 */
-		if (huge+x > one) {
+		if (huge+x > 1.0) {
 			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)));
-	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 */
-		return one + z*(-0.25+(r/s));
+		return 1.0 + z*(-0.25+(r/s));
 	} else {
 		u = 0.5*x;
-		return (one+u)*(one-u) + z*(r/s);
+		return (1.0+u)*(1.0-u) + z*(r/s);
 	}
 }
 
@@ -151,11 +148,11 @@ double y0(double x)
 	ix = 0x7fffffff & hx;
 	/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
 	if (ix >= 0x7ff00000)
-		return one/(x+x*x);
+		return 1.0/(x+x*x);
 	if ((ix|lx) == 0)
-		return -one/zero;
+		return -1.0/0.0;
 	if (hx < 0)
-		return zero/zero;
+		return 0.0/0.0;
 	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
@@ -178,7 +175,7 @@ double y0(double x)
 		 */
 		if (ix < 0x7fe00000) {  /* make sure x+x does not overflow */
 			z = -cos(x+x);
-			if (s*c < zero)
+			if (s*c < 0.0)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -197,7 +194,7 @@ double y0(double x)
 	}
 	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));
 }
 
@@ -286,10 +283,10 @@ static double pzero(double 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;}
-	z = one/(x*x);
+	z = 1.0/(x*x);
 	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;
 }
 
 
@@ -382,8 +379,8 @@ static double qzero(double x)
 	else if (ix >= 0x40122E8B){p = qR5; q = qS5;}
 	else if (ix >= 0x4006DB6D){p = qR3; q = qS3;}
 	else if (ix >= 0x40000000){p = qR2; q = qS2;}
-	z = one/(x*x);
+	z = 1.0/(x*x);
 	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;
 }
diff --git a/src/math/j0f.c b/src/math/j0f.c
index 52cb0ab4..2ee2824b 100644
--- a/src/math/j0f.c
+++ b/src/math/j0f.c
@@ -19,7 +19,6 @@ static float pzerof(float), qzerof(float);
 
 static const float
 huge      = 1e30,
-one       = 1.0,
 invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
 tpi       = 6.3661974669e-01, /* 0x3f22f983 */
 /* R0/S0 on [0, 2.00] */
@@ -32,8 +31,6 @@ S02 =  1.1692678527e-04, /* 0x38f53697 */
 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;
@@ -42,7 +39,7 @@ float j0f(float x)
 	GET_FLOAT_WORD(hx, x);
 	ix = hx & 0x7fffffff;
 	if (ix >= 0x7f800000)
-		return one/(x*x);
+		return 1.0f/(x*x);
 	x = fabsf(x);
 	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
 		s = sinf(x);
@@ -51,7 +48,7 @@ float j0f(float x)
 		cc = s+c;
 		if (ix < 0x7f000000) {  /* make sure x+x does not overflow */
 			z = -cosf(x+x);
-			if (s*c < zero)
+			if (s*c < 0.0f)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -71,20 +68,20 @@ float j0f(float x)
 	}
 	if (ix < 0x39000000) {  /* |x| < 2**-13 */
 		/* raise inexact if x != 0 */
-		if (huge+x > one) {
+		if (huge+x > 1.0f) {
 			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)));
-	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 */
-		return one + z*(-0.25f + (r/s));
+		return 1.0f + z*(-0.25f + (r/s));
 	} else {
 		u = 0.5f*x;
-		return (one+u)*(one-u) + z*(r/s);
+		return (1.0f+u)*(1.0f-u) + z*(r/s);
 	}
 }
 
@@ -110,11 +107,11 @@ float y0f(float x)
 	ix = 0x7fffffff & hx;
 	/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
 	if (ix >= 0x7f800000)
-		return one/(x+x*x);
+		return 1.0f/(x+x*x);
 	if (ix == 0)
-		return -one/zero;
+		return -1.0f/0.0f;
 	if (hx < 0)
-		return zero/zero;
+		return 0.0f/0.0f;
 	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
@@ -137,7 +134,7 @@ float y0f(float x)
 		 */
 		if (ix < 0x7f000000) {  /* make sure x+x not overflow */
 			z = -cosf(x+x);
-			if (s*c < zero)
+			if (s*c < 0.0f)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -156,7 +153,7 @@ float y0f(float x)
 	}
 	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));
 }
 
@@ -244,10 +241,10 @@ static float pzerof(float 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;}
-	z = one/(x*x);
+	z = 1.0f/(x*x);
 	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;
 }
 
 
@@ -340,8 +337,8 @@ static float qzerof(float x)
 	else if (ix >= 0x40f71c58){p = qR5; q = qS5;}
 	else if (ix >= 0x4036db68){p = qR3; q = qS3;}
 	else if (ix >= 0x40000000){p = qR2; q = qS2;}
-	z = one/(x*x);
+	z = 1.0f/(x*x);
 	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;
 }
diff --git a/src/math/j1.c b/src/math/j1.c
index 29ccff0c..4a2cba8d 100644
--- a/src/math/j1.c
+++ b/src/math/j1.c
@@ -60,7 +60,6 @@ static double pone(double), qone(double);
 
 static const double
 huge      = 1e300,
-one       = 1.0,
 invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
 tpi       = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 /* R0/S0 on [0,2] */
@@ -74,8 +73,6 @@ s03 =  1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */
 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;
@@ -84,7 +81,7 @@ double j1(double x)
 	GET_HIGH_WORD(hx, x);
 	ix = hx & 0x7fffffff;
 	if (ix >= 0x7ff00000)
-		return one/x;
+		return 1.0/x;
 	y = fabs(x);
 	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
 		s = sin(y);
@@ -93,7 +90,7 @@ double j1(double x)
 		cc = s-c;
 		if (ix < 0x7fe00000) {  /* make sure y+y not overflow */
 			z = cos(y+y);
-			if (s*c > zero)
+			if (s*c > 0.0)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -116,12 +113,12 @@ double j1(double x)
 	}
 	if (ix < 0x3e400000) {  /* |x| < 2**-27 */
 		/* raise inexact if x!=0 */
-		if (huge+x > one)
+		if (huge+x > 1.0)
 			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;
 }
@@ -151,11 +148,11 @@ double y1(double x)
 	ix = 0x7fffffff & hx;
 	/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 	if (ix >= 0x7ff00000)
-		return one/(x+x*x);
+		return 1.0/(x+x*x);
 	if ((ix|lx) == 0)
-		return -one/zero;
+		return -1.0/0.0;
 	if (hx < 0)
-		return zero/zero;
+		return 0.0/0.0;
 	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
 		s = sin(x);
 		c = cos(x);
@@ -163,7 +160,7 @@ double y1(double x)
 		cc = s-c;
 		if (ix < 0x7fe00000) {  /* make sure x+x not overflow */
 			z = cos(x+x);
-			if (s*c > zero)
+			if (s*c > 0.0)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -192,8 +189,8 @@ double y1(double x)
 		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
@@ -282,10 +279,10 @@ static double pone(double 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;}
-	z = one/(x*x);
+	z = 1.0/(x*x);
 	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
@@ -378,8 +375,8 @@ static double qone(double x)
 	else if (ix >= 0x40122E8B){p = qr5; q = qs5;}
 	else if (ix >= 0x4006DB6D){p = qr3; q = qs3;}
 	else if (ix >= 0x40000000){p = qr2; q = qs2;}
-	z = one/(x*x);
+	z = 1.0/(x*x);
 	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 (.375 + r/s)/x;
 }
diff --git a/src/math/j1f.c b/src/math/j1f.c
index 2d617b67..de459e0d 100644
--- a/src/math/j1f.c
+++ b/src/math/j1f.c
@@ -19,7 +19,6 @@ static float ponef(float), qonef(float);
 
 static const float
 huge      = 1e30,
-one       = 1.0,
 invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */
 tpi       = 6.3661974669e-01, /* 0x3f22f983 */
 /* R0/S0 on [0,2] */
@@ -33,8 +32,6 @@ s03 =  1.1771846857e-06, /* 0x359dffc2 */
 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;
@@ -43,7 +40,7 @@ float j1f(float x)
 	GET_FLOAT_WORD(hx, x);
 	ix = hx & 0x7fffffff;
 	if (ix >= 0x7f800000)
-		return one/x;
+		return 1.0f/x;
 	y = fabsf(x);
 	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
 		s = sinf(y);
@@ -52,7 +49,7 @@ float j1f(float x)
 		cc = s-c;
 		if (ix < 0x7f000000) {  /* make sure y+y not overflow */
 			z = cosf(y+y);
-			if (s*c > zero)
+			if (s*c > 0.0f)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -74,12 +71,12 @@ float j1f(float x)
 	}
 	if (ix < 0x32000000) {  /* |x| < 2**-27 */
 		/* raise inexact if x!=0 */
-		if (huge+x > one)
+		if (huge+x > 1.0f)
 			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;
 }
@@ -108,11 +105,11 @@ float y1f(float x)
 	ix = 0x7fffffff & hx;
 	/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 	if (ix >= 0x7f800000)
-		return one/(x+x*x);
+		return 1.0f/(x+x*x);
 	if (ix == 0)
-		return -one/zero;
+		return -1.0f/0.0f;
 	if (hx < 0)
-		return zero/zero;
+		return 0.0f/0.0f;
 	if (ix >= 0x40000000) {  /* |x| >= 2.0 */
 		s = sinf(x);
 		c = cosf(x);
@@ -120,7 +117,7 @@ float y1f(float x)
 		cc = s-c;
 		if (ix < 0x7f000000) {  /* make sure x+x not overflow */
 			z = cosf(x+x);
-			if (s*c > zero)
+			if (s*c > 0.0f)
 				cc = z/ss;
 			else
 				ss = z/cc;
@@ -149,8 +146,8 @@ float y1f(float x)
 		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
@@ -239,10 +236,10 @@ static float ponef(float 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;}
-	z = one/(x*x);
+	z = 1.0f/(x*x);
 	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
@@ -335,8 +332,8 @@ static float qonef(float x)
 	else if (ix >= 0x40f71c58){p = qr5; q = qs5;}
 	else if (ix >= 0x4036db68){p = qr3; q = qs3;}
 	else if (ix >= 0x40000000){p = qr2; q = qs2;}
-	z = one/(x*x);
+	z = 1.0f/(x*x);
 	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;
 }
diff --git a/src/math/jn.c b/src/math/jn.c
index 082a17bc..d95af15d 100644
--- a/src/math/jn.c
+++ b/src/math/jn.c
@@ -37,12 +37,7 @@
 
 #include "libm.h"
 
-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)
 {
@@ -69,7 +64,7 @@ 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 */
-		b = zero;
+		b = 0.0;
 	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 */
@@ -108,11 +103,11 @@ double jn(int n, double x)
 			 * J(n,x) = 1/n!*(x/2)^n  - ...
 			 */
 			if (n > 33)  /* underflow */
-				b = zero;
+				b = 0.0;
 			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 */
 				}
@@ -165,10 +160,10 @@ double jn(int n, double x)
 				q1 = tmp;
 			}
 			m = n+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);
 			a = t;
-			b = one;
+			b = 1.0;
 			/*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 			 *  Hence, if n*(log(2n/x)) > ...
 			 *  single 8.8722839355e+01
@@ -178,7 +173,7 @@ double jn(int n, double x)
 			 *  likely underflow to zero
 			 */
 			tmp = n;
-			v = two/x;
+			v = 2.0/x;
 			tmp = tmp*log(fabs(v*tmp));
 			if (tmp < 7.09782712893383973096e+02) {
 				for (i=n-1,di=(double)(i+i); i>0; i--) {
@@ -186,7 +181,7 @@ double jn(int n, double x)
 					b *= di;
 					b = b/x - a;
 					a = temp;
-					di -= two;
+					di -= 2.0;
 				}
 			} else {
 				for (i=n-1,di=(double)(i+i); i>0; i--) {
@@ -194,12 +189,12 @@ double jn(int n, double x)
 					b *= di;
 					b = b/x - a;
 					a = temp;
-					di -= two;
+					di -= 2.0;
 					/* scale b to avoid spurious overflow */
 					if (b > 1e100) {
 						a /= b;
 						t /= b;
-						b  = one;
+						b  = 1.0;
 					}
 				}
 			}
@@ -229,9 +224,9 @@ double yn(int n, double x)
 	if ((ix|((uint32_t)(lx|-lx))>>31) > 0x7ff00000)
 		return x+x;
 	if ((ix|lx) == 0)
-		return -one/zero;
+		return -1.0/0.0;
 	if (hx < 0)
-		return zero/zero;
+		return 0.0/0.0;
 	sign = 1;
 	if (n < 0) {
 		n = -n;
@@ -242,7 +237,7 @@ double yn(int n, double x)
 	if (n == 1)
 		return sign*y1(x);
 	if (ix == 0x7ff00000)
-		return zero;
+		return 0.0;
 	if (ix >= 0x52D00000) { /* x > 2**302 */
 		/* (x >> n**2)
 		 *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
diff --git a/src/math/jnf.c b/src/math/jnf.c
index b0b36e6b..fd291103 100644
--- a/src/math/jnf.c
+++ b/src/math/jnf.c
@@ -16,12 +16,6 @@
 #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;
@@ -47,7 +41,7 @@ float jnf(int n, float x)
 	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 */
-		b = zero;
+		b = 0.0f;
 	else if((float)n <= x) {
 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 		a = j0f(x);
@@ -63,11 +57,11 @@ float jnf(int n, float x)
 			 * J(n,x) = 1/n!*(x/2)^n  - ...
 			 */
 			if (n > 33)  /* underflow */
-				b = zero;
+				b = 0.0f;
 			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 */
 				}
@@ -121,10 +115,10 @@ float jnf(int n, float x)
 				q1 = tmp;
 			}
 			m = n+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);
 			a = t;
-			b = one;
+			b = 1.0f;
 			/*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 			 *  Hence, if n*(log(2n/x)) > ...
 			 *  single 8.8722839355e+01
@@ -134,7 +128,7 @@ float jnf(int n, float x)
 			 *  likely underflow to zero
 			 */
 			tmp = n;
-			v = two/x;
+			v = 2.0f/x;
 			tmp = tmp*logf(fabsf(v*tmp));
 			if (tmp < 88.721679688f) {
 				for (i=n-1,di=(float)(i+i); i>0; i--) {
@@ -142,7 +136,7 @@ float jnf(int n, float x)
 					b *= di;
 					b = b/x - a;
 					a = temp;
-					di -= two;
+					di -= 2.0f;
 				}
 			} else {
 				for (i=n-1,di=(float)(i+i); i>0; i--){
@@ -150,12 +144,12 @@ float jnf(int n, float x)
 					b *= di;
 					b = b/x - a;
 					a = temp;
-					di -= two;
+					di -= 2.0f;
 					/* scale b to avoid spurious overflow */
 					if (b > 1e10f) {
 						a /= b;
 						t /= b;
-						b = one;
+						b = 1.0f;
 					}
 				}
 			}
@@ -183,9 +177,9 @@ float ynf(int n, float x)
 	if (ix > 0x7f800000)
 		return x+x;
 	if (ix == 0)
-		return -one/zero;
+		return -1.0f/0.0f;
 	if (hx < 0)
-		return zero/zero;
+		return 0.0f/0.0f;
 	sign = 1;
 	if (n < 0) {
 		n = -n;
@@ -196,7 +190,7 @@ float ynf(int n, float x)
 	if (n == 1)
 		return sign*y1f(x);
 	if (ix == 0x7f800000)
-		return zero;
+		return 0.0f;
 
 	a = y0f(x);
 	b = y1f(x);
diff --git a/src/math/lgamma_r.c b/src/math/lgamma_r.c
index a8ef1956..e3ed1733 100644
--- a/src/math/lgamma_r.c
+++ b/src/math/lgamma_r.c
@@ -82,8 +82,6 @@
 
 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 */
@@ -148,8 +146,6 @@ w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
 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;
@@ -159,7 +155,7 @@ static double sin_pi(double x)
 	ix &= 0x7fffffff;
 
 	if (ix < 0x3fd00000)
-		return __sin(pi*x, zero, 0);
+		return __sin(pi*x, 0.0, 0);
 
 	y = -x;  /* negative x is assumed */
 
@@ -174,7 +170,7 @@ static double sin_pi(double x)
 		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)
@@ -186,14 +182,14 @@ static double sin_pi(double x)
 		}
 	}
 	switch (n) {
-	case 0:  y =  __sin(pi*y, zero, 0); break;
+	case 0:  y =  __sin(pi*y, 0.0, 0); break;
 	case 1:
-	case 2:  y =  __cos(pi*(0.5-y), zero); break;
+	case 2:  y =  __cos(pi*(0.5-y), 0.0); break;
 	case 3:
-	case 4:  y =  __sin(pi*(one-y), zero, 0); break;
+	case 4:  y =  __sin(pi*(1.0-y), 0.0, 0); break;
 	case 5:
-	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;
 	}
 	return -y;
 }
@@ -213,7 +209,7 @@ double __lgamma_r(double x, int *signgamp)
 	if (ix >= 0x7ff00000)
 		return x*x;
 	if ((ix|lx) == 0)
-		return one/zero;
+		return 1.0/0.0;
 	if (ix < 0x3b900000) {  /* |x|<2**-70, return -log(|x|) */
 		if(hx < 0) {
 			*signgamp = -1;
@@ -223,12 +219,12 @@ double __lgamma_r(double x, int *signgamp)
 	}
 	if (hx < 0) {
 		if (ix >= 0x43300000)  /* |x|>=2**52, must be -integer */
-			return one/zero;
+			return 1.0/0.0;
 		t = sin_pi(x);
-		if (t == zero) /* -integer */
-			return one/zero;
+		if (t == 0.0) /* -integer */
+			return 1.0/0.0;
 		nadj = log(pi/fabs(t*x));
-		if (t < zero)
+		if (t < 0.0)
 			*signgamp = -1;
 		x = -x;
 	}
@@ -241,17 +237,17 @@ double __lgamma_r(double x, int *signgamp)
 		if (ix <= 0x3feccccc) {   /* lgamma(x) = lgamma(x+1)-log(x) */
 			r = -log(x);
 			if (ix >= 0x3FE76944) {
-				y = one - x;
+				y = 1.0 - x;
 				i = 0;
 			} else if (ix >= 0x3FCDA661) {
-				y = x - (tc-one);
+				y = x - (tc-1.0);
 				i = 1;
 			} else {
 				y = x;
 				i = 2;
 			}
 		} else {
-			r = zero;
+			r = 0.0;
 			if (ix >= 0x3FFBB4C3) {  /* [1.7316,2] */
 				y = 2.0 - x;
 				i = 0;
@@ -259,7 +255,7 @@ double __lgamma_r(double x, int *signgamp)
 				y = x - tc;
 				i = 1;
 			} else {
-				y = x - one;
+				y = x - 1.0;
 				i = 2;
 			}
 		}
@@ -282,16 +278,16 @@ double __lgamma_r(double x, int *signgamp)
 			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))))));
-		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 */
@@ -303,12 +299,12 @@ double __lgamma_r(double x, int *signgamp)
 		}
 	} else if (ix < 0x43900000) {  /* 8.0 <= x < 2**58 */
 		t = log(x);
-		z = one/x;
+		z = 1.0/x;
 		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.5)*(t-1.0)+w;
 	} else                         /* 2**58 <= x <= inf */
-		r =  x*(log(x)-one);
+		r =  x*(log(x)-1.0);
 	if (hx < 0)
 		r = nadj - r;
 	return r;
diff --git a/src/math/lgammaf_r.c b/src/math/lgammaf_r.c
index f1adcf69..976986aa 100644
--- a/src/math/lgammaf_r.c
+++ b/src/math/lgammaf_r.c
@@ -17,8 +17,6 @@
 
 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 */
@@ -83,8 +81,6 @@ w4  = -5.9518753551e-04, /* 0xba1c065c */
 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;
@@ -109,7 +105,7 @@ static float sin_pif(float x)
 		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)
@@ -125,7 +121,7 @@ static float sin_pif(float x)
 	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;
@@ -148,7 +144,7 @@ float __lgammaf_r(float x, int *signgamp)
 	if (ix >= 0x7f800000)
 		return x*x;
 	if (ix == 0)
-		return one/zero;
+		return 1.0f/0.0f;
 	if (ix < 0x35000000) {  /* |x| < 2**-21, return -log(|x|) */
 		if (hx < 0) {
 			*signgamp = -1;
@@ -158,12 +154,12 @@ float __lgammaf_r(float x, int *signgamp)
 	}
 	if (hx < 0) {
 		if (ix >= 0x4b000000)  /* |x| >= 2**23, must be -integer */
-			return one/zero;
+			return 1.0f/0.0f;
 		t = sin_pif(x);
-		if (t == zero) /* -integer */
-			return one/zero;
+		if (t == 0.0f) /* -integer */
+			return 1.0f/0.0f;
 		nadj = logf(pi/fabsf(t*x));
-		if (t < zero)
+		if (t < 0.0f)
 			*signgamp = -1;
 		x = -x;
 	}
@@ -176,17 +172,17 @@ float __lgammaf_r(float x, int *signgamp)
 		if (ix <= 0x3f666666) {  /* lgamma(x) = lgamma(x+1)-log(x) */
 			r = -logf(x);
 			if (ix >= 0x3f3b4a20) {
-				y = one - x;
+				y = 1.0f - x;
 				i = 0;
 			} else if (ix >= 0x3e6d3308) {
-				y = x - (tc-one);
+				y = x - (tc-1.0f);
 				i = 1;
 			} else {
 				y = x;
 				i = 2;
 			}
 		} else {
-			r = zero;
+			r = 0.0f;
 			if (ix >= 0x3fdda618) {  /* [1.7316,2] */
 				y = 2.0f - x;
 				i = 0;
@@ -194,7 +190,7 @@ float __lgammaf_r(float x, int *signgamp)
 				y = x - tc;
 				i = 1;
 			} else {
-				y = x - one;
+				y = x - 1.0f;
 				i = 2;
 			}
 		}
@@ -217,16 +213,16 @@ float __lgammaf_r(float x, int *signgamp)
 			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))))));
-		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 */
@@ -238,12 +234,12 @@ float __lgammaf_r(float x, int *signgamp)
 		}
 	} else if (ix < 0x5c800000) {  /* 8.0 <= x < 2**58 */
 		t = logf(x);
-		z = one/x;
+		z = 1.0f/x;
 		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 */
-		r =  x*(logf(x)-one);
+		r =  x*(logf(x)-1.0f);
 	if (hx < 0)
 		r = nadj - r;
 	return r;
diff --git a/src/math/lgammal.c b/src/math/lgammal.c
index ec7c9a04..8fae1be8 100644
--- a/src/math/lgammal.c
+++ b/src/math/lgammal.c
@@ -95,8 +95,6 @@ long double __lgammal_r(long double x, int *sg)
 }
 #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,
 
@@ -200,8 +198,6 @@ w5 =  8.412723297322498080632E-4L,
 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;
@@ -226,7 +222,7 @@ static long double sin_pi(long double x)
 		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 */
@@ -244,11 +240,11 @@ static long double sin_pi(long double x)
 		break;
 	case 1:
 	case 2:
-		y = cosl(pi * (half - y));
+		y = cosl(pi * (0.5 - y));
 		break;
 	case 3:
 	case 4:
-		y = sinl(pi * (one - y));
+		y = sinl(pi * (1.0 - y));
 		break;
 	case 5:
 	case 6:
@@ -273,7 +269,7 @@ long double __lgammal_r(long double x, int *sg) {
 	if ((ix | i0 | i1) == 0) {
 		if (se & 0x8000)
 			*sg = -1;
-		return one / fabsl(x);
+		return 1.0 / fabsl(x);
 	}
 
 	ix = (ix << 16) | (i0 >> 16);
@@ -291,10 +287,10 @@ long double __lgammal_r(long double x, int *sg) {
 	}
 	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));
-		if (t < zero)
+		if (t < 0.0)
 			*sg = -1;
 		x = -x;
 	}
@@ -306,19 +302,19 @@ long double __lgammal_r(long double x, int *sg) {
 	else if (ix < 0x40008000) {  /* x < 2.0 */
 		if (ix <= 0x3ffee666) {  /* 8.99993896484375e-1 */
 			/* lgamma(x) = lgamma(x+1) - log(x) */
-			r = -logl (x);
+			r = -logl(x);
 			if (ix >= 0x3ffebb4a) {  /* 7.31597900390625e-1 */
-				y = x - one;
+				y = x - 1.0;
 				i = 0;
 			} else if (ix >= 0x3ffced33) {  /* 2.31639862060546875e-1 */
-				y = x - (tc - one);
+				y = x - (tc - 1.0);
 				i = 1;
 			} else { /* x < 0.23 */
 				y = x;
 				i = 2;
 			}
 		} else {
-			r = zero;
+			r = 0.0;
 			if (ix >= 0x3fffdda6) {  /* 1.73162841796875 */
 				/* [1.7316,2] */
 				y = x - 2.0;
@@ -329,7 +325,7 @@ long double __lgammal_r(long double x, int *sg) {
 				i = 1;
 			} else {
 				/* [0.9, 1.23] */
-				y = x - one;
+				y = x - 1.0;
 				i = 2;
 			}
 		}
@@ -337,7 +333,7 @@ long double __lgammal_r(long double x, int *sg) {
 		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)))));
@@ -348,17 +344,17 @@ long double __lgammal_r(long double x, int *sg) {
 		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;
-		t = zero;
+		t = 0.0;
 		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 */
@@ -370,18 +366,18 @@ long double __lgammal_r(long double x, int *sg) {
 			z *= (y + 3.0); /* FALLTHRU */
 		case 3:
 			z *= (y + 2.0); /* FALLTHRU */
-			r += logl (z);
+			r += logl(z);
 			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))))));
-		r = (x - half) * (t - one) + w;
+		r = (x - 0.5) * (t - 1.0) + w;
 	} else /* 2**66 <= x <= inf */
-		r = x * (logl (x) - one);
+		r = x * (logl(x) - 1.0);
 	if (se & 0x8000)
 		r = nadj - r;
 	return r;
diff --git a/src/math/log.c b/src/math/log.c
index 1bb006a3..98051205 100644
--- a/src/math/log.c
+++ b/src/math/log.c
@@ -74,8 +74,6 @@ Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
 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;
@@ -87,9 +85,9 @@ double log(double x)
 	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 */
 		if (hx < 0)
-			return (x-x)/zero;   /* log(-#) = NaN */
+			return (x-x)/0.0;   /* log(-#) = NaN */
 		/* subnormal number, scale up x */
 		k -= 54;
 		x *= two54;
@@ -104,9 +102,9 @@ double log(double x)
 	k += i>>20;
 	f = x - 1.0;
 	if ((0x000fffff&(2+hx)) < 3) {  /* -2**-20 <= f < 2**-20 */
-		if (f == zero) {
+		if (f == 0.0) {
 			if (k == 0) {
-				return zero;
+				return 0.0;
 			}
 			dk = (double)k;
 			return dk*ln2_hi + dk*ln2_lo;
diff --git a/src/math/log10.c b/src/math/log10.c
index 5422599a..ed65d9be 100644
--- a/src/math/log10.c
+++ b/src/math/log10.c
@@ -27,8 +27,6 @@ ivln10lo  = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
 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;
@@ -40,9 +38,9 @@ double log10(double x)
 	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 */
 		if (hx<0)
-			return (x-x)/zero;   /* log(-#) = NaN */
+			return (x-x)/0.0;   /* log(-#) = NaN */
 		/* subnormal number, scale up x */
 		k -= 54;
 		x *= two54;
@@ -51,7 +49,7 @@ double log10(double x)
 	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;
diff --git a/src/math/log10f.c b/src/math/log10f.c
index 7e4ac9a8..e10749b5 100644
--- a/src/math/log10f.c
+++ b/src/math/log10f.c
@@ -23,8 +23,6 @@ ivln10lo  = -3.1689971365e-05, /* 0xb804ead9 */
 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;
@@ -35,9 +33,9 @@ float log10f(float x)
 	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 */
 		if (hx < 0)
-			return (x-x)/zero;   /* log(-#) = NaN */
+			return (x-x)/0.0f;   /* log(-#) = NaN */
 		/* subnormal number, scale up x */
 		k -= 25;
 		x *= two25;
@@ -46,7 +44,7 @@ float log10f(float x)
 	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;
diff --git a/src/math/log1p.c b/src/math/log1p.c
index f7154d0c..6c67249c 100644
--- a/src/math/log1p.c
+++ b/src/math/log1p.c
@@ -88,8 +88,6 @@ Lp5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
 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;
@@ -102,12 +100,12 @@ double log1p(double x)
 	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 */
-			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;
 		}
@@ -151,9 +149,9 @@ double log1p(double x)
 	}
 	hfsq = 0.5*f*f;
 	if (hu == 0) {   /* |f| < 2**-20 */
-		if (f == zero) {
+		if (f == 0.0) {
 			if(k == 0)
-				return zero;
+				return 0.0;
 			c += k*ln2_lo;
 			return k*ln2_hi + c;
 		}
diff --git a/src/math/log1pf.c b/src/math/log1pf.c
index 75eeb371..39832d28 100644
--- a/src/math/log1pf.c
+++ b/src/math/log1pf.c
@@ -27,8 +27,6 @@ Lp5 = 1.8183572590e-01, /* 3E3A3325 */
 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;
@@ -41,12 +39,12 @@ float log1pf(float x)
 	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 */
-			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;
 		}
@@ -91,9 +89,9 @@ float log1pf(float x)
 	}
 	hfsq = 0.5f * f * f;
 	if (hu == 0) {  /* |f| < 2**-20 */
-		if (f == zero) {
+		if (f == 0.0f) {
 			if (k == 0)
-				return zero;
+				return 0.0f;
 			c += k*ln2_lo;
 			return k*ln2_hi+c;
 		}
diff --git a/src/math/log1pl.c b/src/math/log1pl.c
index 17eb4cef..1400d365 100644
--- a/src/math/log1pl.c
+++ b/src/math/log1pl.c
@@ -118,11 +118,11 @@ long double log1pl(long double xm1)
 	if (xm1 == 0.0)
 		return xm1;
 
-	x = xm1 + 1.0L;
+	x = xm1 + 1.0;
 
 	/* Test for domain errors.  */
-	if (x <= 0.0L) {
-		if (x == 0.0L)
+	if (x <= 0.0) {
+		if (x == 0.0)
 			return -INFINITY;
 		return NAN;
 	}
@@ -136,12 +136,12 @@ long double log1pl(long double xm1)
 	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)   */
-			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;
 		}
 		x = z / y;
 		z = x*x;
@@ -156,12 +156,12 @@ long double log1pl(long double xm1)
 	if (x < SQRTH) {
 		e -= 1;
 		if (e != 0)
-			x = 2.0 * x - 1.0L;
+			x = 2.0 * x - 1.0;
 		else
 			x = xm1;
 	} else {
 		if (e != 0)
-			x = x - 1.0L;
+			x = x - 1.0;
 		else
 			x = xm1;
 	}
diff --git a/src/math/log2.c b/src/math/log2.c
index a5b8abdd..1974215d 100644
--- a/src/math/log2.c
+++ b/src/math/log2.c
@@ -27,8 +27,6 @@ two54   = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
 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;
@@ -40,9 +38,9 @@ double log2(double x)
 	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 */
 		if (hx < 0)
-			return (x-x)/zero;   /* log(-#) = NaN */
+			return (x-x)/0.0;   /* log(-#) = NaN */
 		/* subnormal number, scale up x */
 		k -= 54;
 		x *= two54;
@@ -51,7 +49,7 @@ double log2(double x)
 	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;
diff --git a/src/math/log2f.c b/src/math/log2f.c
index c9b9282c..e0d6a9e4 100644
--- a/src/math/log2f.c
+++ b/src/math/log2f.c
@@ -21,8 +21,6 @@ two25   =  3.3554432000e+07, /* 0x4c000000 */
 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;
@@ -33,9 +31,9 @@ float log2f(float x)
 	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 */
 		if (hx < 0)
-			return (x-x)/zero;   /* log(-#) = NaN */
+			return (x-x)/0.0f;   /* log(-#) = NaN */
 		/* subnormal number, scale up x */
 		k -= 25;
 		x *= two25;
@@ -44,7 +42,7 @@ float log2f(float x)
 	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;
diff --git a/src/math/logf.c b/src/math/logf.c
index a4ed697b..c7f7dbe6 100644
--- a/src/math/logf.c
+++ b/src/math/logf.c
@@ -25,8 +25,6 @@ Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
 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;
@@ -37,9 +35,9 @@ float logf(float x)
 	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 */
 		if (ix < 0)
-			return (x-x)/zero;   /* log(-#) = NaN */
+			return (x-x)/0.0f;   /* log(-#) = NaN */
 		/* subnormal number, scale up x */
 		k -= 25;
 		x *= two25;
@@ -54,9 +52,9 @@ float logf(float x)
 	k += i>>23;
 	f = x - 1.0f;
 	if ((0x007fffff & (0x8000 + ix)) < 0xc000) {  /* -2**-9 <= f < 2**-9 */
-		if (f == zero) {
+		if (f == 0.0f) {
 			if (k == 0)
-				return zero;
+				return 0.0f;
 			dk = (float)k;
 			return dk*ln2_hi + dk*ln2_lo;
 		}
diff --git a/src/math/modf.c b/src/math/modf.c
index ff85b2a3..cca3b652 100644
--- a/src/math/modf.c
+++ b/src/math/modf.c
@@ -21,8 +21,6 @@
 
 #include "libm.h"
 
-static const double one = 1.0;
-
 double modf(double x, double *iptr)
 {
 	int32_t i0,i1,j0;
@@ -51,7 +49,7 @@ double modf(double x, double *iptr)
 			*iptr = x;
 			return 0.0 / x;
 		}
-		*iptr = x*one;
+		*iptr = x;
 		GET_HIGH_WORD(high, x);
 		INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
 		return x;
diff --git a/src/math/pow.c b/src/math/pow.c
index 5bcedd5b..052825a6 100644
--- a/src/math/pow.c
+++ b/src/math/pow.c
@@ -59,9 +59,6 @@ static const double
 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,
@@ -101,15 +98,15 @@ double pow(double x, double y)
 	ix = hx & 0x7fffffff;
 	iy = hy & 0x7fffffff;
 
-	/* y == zero: x**0 = 1 */
+	/* y == 0.0: x**0 = 1 */
 	if ((iy|ly) == 0)
-		return one;
+		return 1.0;
 
 	/* x == 1: 1**y = 1, even if y is NaN */
 	if (hx == 0x3ff00000 && lx == 0)
-		return one;
+		return 1.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);
@@ -141,15 +138,15 @@ double pow(double x, double y)
 	if (ly == 0) {
 		if (iy == 0x7ff00000) {  /* y is +-inf */
 			if (((ix-0x3ff00000)|lx) == 0)  /* (-1)**+-inf is 1 */
-				return one;
+				return 1.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 */
-				return hy < 0 ? -y : zero;
+				return hy < 0 ? -y : 0.0;
 		}
 		if (iy == 0x3ff00000) {  /* y is +-1 */
 			if (hy < 0)
-				return one/x;
+				return 1.0/x;
 			return x;
 		}
 		if (hy == 0x40000000)    /* y is 2 */
@@ -166,7 +163,7 @@ double pow(double x, double y)
 		if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */
 			z = ax;
 			if (hy < 0)   /* z = (1/|x|) */
-				z = one/z;
+				z = 1.0/z;
 			if (hx < 0) {
 				if (((ix-0x3ff00000)|yisint) == 0) {
 					z = (z-z)/(z-z); /* (-1)**non-int is NaN */
@@ -187,9 +184,9 @@ double pow(double x, double y)
 	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 */
 	if ((n|(yisint-1)) == 0)
-		s = -one;/* (-ve)**(odd int) */
+		s = -1.0;/* (-ve)**(odd int) */
 
 	/* |y| is huge */
 	if (iy > 0x41e00000) { /* if |y| > 2**31 */
@@ -206,7 +203,7 @@ double pow(double x, double y)
 			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;
@@ -239,12 +236,12 @@ double pow(double x, double y)
 
 		/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
 		u = ax - bp[k];        /* bp[0]=1.0, bp[1]=1.5 */
-		v = one/(ax+bp[k]);
+		v = 1.0/(ax+bp[k]);
 		ss = u*v;
 		s_h = ss;
 		SET_LOW_WORD(s_h, 0);
 		/* t_h=ax+bp[k] High */
-		t_h = zero;
+		t_h = 0.0;
 		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);
@@ -299,7 +296,7 @@ double pow(double x, double y)
 	if (i > 0x3fe00000) {  /* if |z| > 0.5, set n = [z+0.5] */
 		n = j + (0x00100000>>(k+1));
 		k = ((n&0x7fffffff)>>20) - 0x3ff;  /* new k for n */
-		t = zero;
+		t = 0.0;
 		SET_HIGH_WORD(t, n & ~(0x000fffff>>k));
 		n = ((n&0x000fffff)|0x00100000)>>(20-k);
 		if (j < 0)
@@ -314,8 +311,8 @@ double pow(double x, double y)
 	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 */
diff --git a/src/math/powf.c b/src/math/powf.c
index 01aced0e..c101d95c 100644
--- a/src/math/powf.c
+++ b/src/math/powf.c
@@ -19,9 +19,6 @@ static const float
 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,
@@ -60,15 +57,15 @@ float powf(float x, float y)
 	ix = hx & 0x7fffffff;
 	iy = hy & 0x7fffffff;
 
-	/* y == zero: x**0 = 1 */
+	/* y == 0: x**0 = 1 */
 	if (iy == 0)
-		return one;
+		return 1.0f;
 
 	/* x == 1: 1**y = 1, even if y is NaN */
 	if (hx == 0x3f800000)
-		return one;
+		return 1.0f;
 
-	/* 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);
 
@@ -92,15 +89,15 @@ float powf(float x, float y)
 	/* special value of y */
 	if (iy == 0x7f800000) {  /* y is +-inf */
 		if (ix == 0x3f800000)      /* (-1)**+-inf is 1 */
-			return one;
+			return 1.0f;
 		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 */
-			return hy < 0 ? -y : zero;
+			return hy < 0 ? -y : 0.0f;
 	}
 	if (iy == 0x3f800000) {  /* y is +-1 */
 		if (hy < 0)
-			return one/x;
+			return 1.0f/x;
 		return x;
 	}
 	if (hy == 0x40000000)    /* y is 2 */
@@ -115,7 +112,7 @@ float powf(float x, float y)
 	if (ix == 0x7f800000 || ix == 0 || ix == 0x3f800000) { /* x is +-0,+-inf,+-1 */
 		z = ax;
 		if (hy < 0)  /* z = (1/|x|) */
-			z = one/z;
+			z = 1.0f/z;
 		if (hx < 0) {
 			if (((ix-0x3f800000)|yisint) == 0) {
 				z = (z-z)/(z-z); /* (-1)**non-int is NaN */
@@ -131,9 +128,9 @@ float powf(float x, float y)
 	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) */
-		sn = -one;
+		sn = -1.0f;
 
 	/* |y| is huge */
 	if (iy > 0x4d000000) { /* if |y| > 2**27 */
@@ -178,7 +175,7 @@ float powf(float x, float y)
 
 		/* 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 */
-		v = one/(ax+bp[k]);
+		v = 1.0f/(ax+bp[k]);
 		s = u*v;
 		s_h = s;
 		GET_FLOAT_WORD(is, s_h);
@@ -257,8 +254,8 @@ float powf(float x, float y)
 	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 */
diff --git a/src/math/remquol.c b/src/math/remquol.c
index dd18f35c..721231b4 100644
--- a/src/math/remquol.c
+++ b/src/math/remquol.c
@@ -48,7 +48,7 @@ typedef uint32_t manh_t;
 
 #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
diff --git a/src/math/rintl.c b/src/math/rintl.c
index 1cc35df5..6a311a9d 100644
--- a/src/math/rintl.c
+++ b/src/math/rintl.c
@@ -79,7 +79,7 @@ long double rintl(long double x)
 	 * 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;
diff --git a/src/math/scalbnl.c b/src/math/scalbnl.c
index a5d0adba..c605b8da 100644
--- a/src/math/scalbnl.c
+++ b/src/math/scalbnl.c
@@ -29,7 +29,7 @@ long double scalbnl(long double x, int n)
 				return x * 0x1p-16382L;
 		}
 	}
-	scale.e = 1.0L;
+	scale.e = 1.0;
 	scale.bits.exp = 0x3fff + n;
 	return x * scale.e;
 }
diff --git a/src/math/sinh.c b/src/math/sinh.c
index 935879c5..0c67ba39 100644
--- a/src/math/sinh.c
+++ b/src/math/sinh.c
@@ -29,7 +29,7 @@
 
 #include "libm.h"
 
-static const double one = 1.0, huge = 1.0e307;
+static const double huge = 1.0e307;
 
 double sinh(double x)
 {
@@ -50,12 +50,12 @@ double sinh(double x)
 	if (ix < 0x40360000) {  /* |x|<22 */
 		if (ix < 0x3e300000)  /* |x|<2**-28 */
 			/* raise inexact, return x */
-			if (huge+x > one)
+			if (huge+x > 1.0)
 				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|) */
diff --git a/src/math/sinhf.c b/src/math/sinhf.c
index fd11b849..b8d88224 100644
--- a/src/math/sinhf.c
+++ b/src/math/sinhf.c
@@ -15,7 +15,7 @@
 
 #include "libm.h"
 
-static const float one = 1.0, huge = 1.0e37;
+static const float huge = 1.0e37;
 
 float sinhf(float x)
 {
@@ -36,12 +36,12 @@ float sinhf(float x)
 	if (ix < 0x41100000) {   /* |x|<9 */
 		if (ix < 0x39800000)  /* |x|<2**-12 */
 			/* raise inexact, return x */
-			if (huge+x > one)
+			if (huge+x > 1.0f)
 				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|) */
diff --git a/src/math/sinhl.c b/src/math/sinhl.c
index 2252dec9..8a54677e 100644
--- a/src/math/sinhl.c
+++ b/src/math/sinhl.c
@@ -35,7 +35,7 @@ long double sinhl(long double x)
 	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)
 {
@@ -55,12 +55,12 @@ 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 */
-			if (huge + x > one)
+			if (huge + x > 1.0)
 				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|) */
diff --git a/src/math/sqrt.c b/src/math/sqrt.c
index 2ebd022b..b2775673 100644
--- a/src/math/sqrt.c
+++ b/src/math/sqrt.c
@@ -78,7 +78,7 @@
 
 #include "libm.h"
 
-static const double one = 1.0, tiny = 1.0e-300;
+static const double tiny = 1.0e-300;
 
 double sqrt(double x)
 {
@@ -161,13 +161,13 @@ double sqrt(double x)
 
 	/* 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++;
-			} else if (z > one) {
+			} else if (z > 1.0) {
 				if (q1 == (uint32_t)0xfffffffe)
 					q++;
 				q1 += 2;
diff --git a/src/math/sqrtf.c b/src/math/sqrtf.c
index 35c24e50..28cb4ad3 100644
--- a/src/math/sqrtf.c
+++ b/src/math/sqrtf.c
@@ -15,7 +15,7 @@
 
 #include "libm.h"
 
-static const float one = 1.0, tiny = 1.0e-30;
+static const float tiny = 1.0e-30;
 
 float sqrtf(float x)
 {
@@ -68,10 +68,10 @@ float sqrtf(float x)
 
 	/* 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)
 				q += 2;
 			else
 				q += q & 1;
diff --git a/src/math/tanh.c b/src/math/tanh.c
index 957c43e9..21138643 100644
--- a/src/math/tanh.c
+++ b/src/math/tanh.c
@@ -35,7 +35,7 @@
 
 #include "libm.h"
 
-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;
 
 double tanh(double x)
 {
@@ -48,26 +48,26 @@ double tanh(double x)
 	/* 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 */
 		else
-			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 (huge+x > one)
+			if (huge+x > 1.0f)
 				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);
 		} else {
-			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 */
-		z = one - tiny;  /* raise inexact */
+		z = 1.0f - tiny;  /* raise inexact */
 	}
 	return jx >= 0 ? z : -z;
 }
diff --git a/src/math/tanhf.c b/src/math/tanhf.c
index 97d0eb53..7cb459d0 100644
--- a/src/math/tanhf.c
+++ b/src/math/tanhf.c
@@ -15,7 +15,9 @@
 
 #include "libm.h"
 
-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;
 
 float tanhf(float x)
 {
@@ -28,26 +30,26 @@ float tanhf(float x)
 	/* 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 */
 		else
-			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 (huge+x > one)
+			if (huge+x > 1.0f)
 				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);
 		} else {
-			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 */
-		z = one - tiny;  /* raise inexact */
+		z = 1.0f - tiny;  /* raise inexact */
 	}
 	return jx >= 0 ? z : -z;
 }
diff --git a/src/math/tanhl.c b/src/math/tanhl.c
index e62be59b..92efb20d 100644
--- a/src/math/tanhl.c
+++ b/src/math/tanhl.c
@@ -41,7 +41,7 @@ long double tanhl(long double x)
 	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)
 {
@@ -57,8 +57,8 @@ long double tanhl(long double x)
 	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 */
 	}
 
 	/* |x| < 23 */
@@ -66,17 +66,17 @@ long double tanhl(long double x)
 		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  */
-			t = expm1l(two*fabsl(x));
-			z = one - two/(t+two);
+			t = expm1l(2.0*fabsl(x));
+			z = 1.0 - 2.0/(t+2.0);
 		} else {
-			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 {
-		z = one - tiny;  /* raise inexact flag */
+		z = 1.0 - tiny;  /* raise inexact flag */
 	}
 	return se & 0x8000 ? -z : z;
 }
diff --git a/src/math/tgammal.c b/src/math/tgammal.c
index 1b8fddea..1b460da5 100644
--- a/src/math/tgammal.c
+++ b/src/math/tgammal.c
@@ -183,18 +183,18 @@ static long double stirf(long double x)
 {
 	long double y, w, v;
 
-	w = 1.0L/x;
+	w = 1.0/x;
 	/* For large x, use rational coefficients from the analytical expansion.  */
-	if (x > 1024.0L)
+	if (x > 1024.0)
 		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
-		 + 1.0L;
+		 + 1.0;
 	else
-		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);
@@ -219,10 +219,10 @@ long double tgammal(long double x)
 	if (x == -INFINITY)
 		return x - x;
 	q = fabsl(x);
-	if (q > 13.0L) {
+	if (q > 13.0) {
 		if (q > MAXGAML)
 			goto goverf;
-		if (x < 0.0L) {
+		if (x < 0.0) {
 			p = floorl(q);
 			if (p == q)
 				return (x - x) / (x - x);
@@ -231,7 +231,7 @@ long double tgammal(long double x)
 				signgam = -1;
 			z = q - p;
 			if (z > 0.5L) {
-				p += 1.0L;
+				p += 1.0;
 				z = q - p;
 			}
 			z = q * sinl(PIL * z);
@@ -247,25 +247,25 @@ goverf:
 		return signgam * 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;
-		x += 1.0L;
+		x += 1.0;
 	}
 	if (x <= 0.03125L)
 		goto small;
-	while (x < 2.0L) {
+	while (x < 2.0) {
 		z /= x;
-		x += 1.0L;
+		x += 1.0;
 	}
-	if (x == 2.0L)
+	if (x == 2.0)
 		return z;
 
-	x -= 2.0L;
+	x -= 2.0;
 	p = __polevll(x, P, 7);
 	q = __polevll(x, Q, 8);
 	z = z * p / q;
@@ -274,9 +274,9 @@ goverf:
 	return z;
 
 small:
-	if (x == 0.0L)
+	if (x == 0.0)
 		return (x - x) / (x - x);
-	if (x < 0.0L) {
+	if (x < 0.0) {
 		x = -x;
 		q = z / (x * __polevll(x, SN, 8));
 		signgam = -1;
-- 
cgit v1.2.1