svn commit: r240864 - in head/lib/msun: ld128 ld80

[Steve Kargl]

Author: kargl
Date: Sun Sep 23 18:06:27 2012
New Revision: 240864
URL: http://svn.freebsd.org/changeset/base/240864

Log:
  * ld80/s_expl.c:
    . Guard a comment from reformatting by indent(1).
    . Re-order variables in declarations to alphabetical order.
    . Remove a banal comment.
  
  * ld128/s_expl.c:
    . Add a comment to point to ld80/s_expl.c for implementation details.
    . Move the #define of INTERVAL to reduce the diff with ld80/s_expl.c.
    . twom10000 does not need to be volatile, so move its declaration.
    . Re-order variables in declarations to alphabetical order.
    . Add a comment that describes the argument reduction.
    . Remove the same banal comment found in ld80/s_expl.c.
  
  Reviewed by:	bde
  Approved by:	das (mentor)

Modified:
  head/lib/msun/ld128/s_expl.c
  head/lib/msun/ld80/s_expl.c

Modified: head/lib/msun/ld128/s_expl.c
==============================================================================
--- head/lib/msun/ld128/s_expl.c	Sun Sep 23 17:47:01 2012	(r240863)
+++ head/lib/msun/ld128/s_expl.c	Sun Sep 23 18:06:27 2012	(r240864)
@@ -27,24 +27,29 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * ld128 version of s_expl.c.  See ../ld80/s_expl.c for most comments.
+ */
+
 #include <float.h>
 
 #include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
 
+#define	INTERVALS	128
 #define	BIAS	(LDBL_MAX_EXP - 1)
 
-/* XXX Prevent gcc from erroneously constant folding this: */
-static volatile const long double twom10000 = 0x1p-10000L, tiny = 0x1p-10000L;
+static volatile const long double tiny = 0x1p-10000L;
 
 static const long double
-huge = 0x1p10000L,
-o_threshold =  11356.523406294143949491931077970763428L,
-u_threshold = -11433.462743336297878837243843452621503L,
+INV_L = 1.84664965233787316142070359168242182e+02L,
 L1 = 5.41521234812457272982212595914567508e-03L,
 L2 = -1.02536706388947310094527932552595546e-29L,
-INV_L = 1.84664965233787316142070359168242182e+02L;
+huge = 0x1p10000L,
+o_threshold =  11356.523406294143949491931077970763428L,
+twom10000 = 0x1p-10000L,
+u_threshold = -11433.462743336297878837243843452621503L;
 
 static const long double
 P2 = 5.00000000000000000000000000000000000e-1L,
@@ -58,8 +63,6 @@ P9 = 2.755731922401038678178761995444688
 P10 = 2.75573236172670046201884000197885520e-7L,
 P11 = 2.50517544183909126492878226167697856e-8L;
 
-#define	INTERVALS	128
-
 static const struct {
 	long double	hi;
 	long double	lo;
@@ -225,9 +228,10 @@ expl(long double x)
 			return (1.0L + x);
 	}
 
+	/* Reduce x to (k*ln2 + midpoint[n2] + r1 + r2). */
 	fn = x * INV_L + 0x1.8p112 - 0x1.8p112;
 	n  = (int)fn;
-	n2 = (unsigned)n % INTERVALS;		/* Tang's j. */
+	n2 = (unsigned)n % INTERVALS;
 	k = (n - n2) / INTERVALS;
 	r1 = x - fn * L1;
 	r2 = -fn * L2;

Modified: head/lib/msun/ld80/s_expl.c
==============================================================================
--- head/lib/msun/ld80/s_expl.c	Sun Sep 23 17:47:01 2012	(r240863)
+++ head/lib/msun/ld80/s_expl.c	Sun Sep 23 18:06:27 2012	(r240864)
@@ -29,7 +29,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/*
+/*-
  * Compute the exponential of x for Intel 80-bit format.  This is based on:
  *
  *   PTP Tang, "Table-driven implementation of the exponential function
@@ -70,9 +70,9 @@ static const double __aligned(64)
  * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
  * bits zero so that multiplication of it by n is exact.
  */
+INV_L = 1.8466496523378731e+2,		/*  0x171547652b82fe.0p-45 */
 L1 =  5.4152123484527692e-3,		/*  0x162e42ff000000.0p-60 */
 L2 = -3.2819649005320973e-13,		/* -0x1718432a1b0e26.0p-94 */
-INV_L = 1.8466496523378731e+2,		/*  0x171547652b82fe.0p-45 */
 /*
  * Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]:
  * |exp(x) - p(x)| < 2**-77.2
@@ -231,7 +231,7 @@ long double
 expl(long double x)
 {
 	union IEEEl2bits u, v;
-	long double fn, r, r1, r2, q, t, t23, t45, twopk, twopkp10000, z;
+	long double fn, q, r, r1, r2, t, t23, t45, twopk, twopkp10000, z;
 	int k, n, n2;
 	uint16_t hx, ix;
 
@@ -268,7 +268,7 @@ expl(long double x)
 #else
 	n  = (int)fn;
 #endif
-	n2 = (unsigned)n % INTERVALS;		/* Tang's j. */
+	n2 = (unsigned)n % INTERVALS;
 	k = (n - n2) / INTERVALS;
 	r1 = x - fn * L1;
 	r2 = -fn * L2;