Merge remote branch 'nsz/master'

94 files changed, 654 insertions, 806 deletions

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.h b/src/math/__rem_pio2l.c
index 11123c3f..f0cb99ac 100644
--- a/src/math/__rem_pio2l.h
+++ 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 */
@@ -44,7 +43,7 @@ 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)
+int __rem_pio2l(long double x, long double *y)
 {
 	union IEEEl2bits u,u1;
 	long double z,w,t,r,fn;
@@ -114,7 +113,7 @@ static inline 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/__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;
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/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,
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/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/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/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/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/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/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
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/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/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/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/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/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/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/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. */
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/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;
 }
 
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) {
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/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/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/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/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/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;
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;