Implementation of expl()

David Schultz das at FreeBSD.ORG
Sun Dec 9 14:42:01 PST 2007


On Sun, Dec 09, 2007, David Schultz wrote:
> On Sun, Dec 09, 2007, Steve Kargl wrote:
> > On Sun, Dec 09, 2007 at 04:25:05PM -0500, David Schultz wrote:
> > > On Fri, Nov 02, 2007, Steve Kargl wrote:
> > > > With all the success of developing several other missing
> > > > C99 math routines, I decided to tackle expl() (which I
> > > > need to tackle tanhl()).
> > > 
> > > Hmm, great, but where's the patch? :)  Maybe the mailing list
> > > software ate it.
> > 
> > This is the current version.  I need to revise how I computed
> > the ploynomial coefficient.  In short, I mapped r in 
> > [-ln(2)/2:ln(2)/2] into the range x in [-1,1] for the Chebyshev
> > interpolation, but I never scaled x back into r.  This is the
> > reason why the lines "r = r * TOLN2;" exists.
> 
> Some minor fixes (mostly whitespace) are below.

Oops, here it is.

--- /usr/src/lib/msun/src/e_expl.c.bak	2007-12-09 17:30:35.000000000 -0500
+++ /usr/src/lib/msun/src/e_expl.c	2007-12-09 17:34:14.000000000 -0500
@@ -38,45 +38,46 @@
 #define XMAX	0x2.c5c85fdf473de6ap12L
 
 /* ln(LDBL_MIN) = -11355.137111933024 */
-#define XMIN   -0x2.c5b2319c4843accp12L
+#define XMIN	-0x2.c5b2319c4843accp12L
 
 /* | ln(smallest subnormal) | = 11399.498531488861 */
 #define GRAD	0x2.c877f9fc278aeaap12L
 
-#define LN2      0xb.17217f7d1cf79acp-4L /* ln(2) */
-#define LN2HI    0xb.17217f7d2000000p-4L
-#define LN2LO   -0x3.08654361c4c67fcp-44L
+#define LN2	0xb.17217f7d1cf79acp-4L /* ln(2) */
+#define LN2HI	0xb.17217f7d2000000p-4L
+#define LN2LO	-0x3.08654361c4c67fcp-44L
 
-#define TOLN2    0x2.e2a8eca5705fc2fp0L  /* 2/ln(2) */
+#define TOLN2	0x2.e2a8eca5705fc2fp0L  /* 2/ln(2) */
 
-#define ILN2     0x1.71547652b82fe17p0L  /* 1/ln(2) */
-#define ILN2HI   0x1.71547652b800000p0L
-#define ILN2LO   0x2.fe1777d0ffda0d2p-44L
+#define ILN2	0x1.71547652b82fe17p0L  /* 1/ln(2) */
+#define ILN2HI	0x1.71547652b800000p0L
+#define ILN2LO	0x2.fe1777d0ffda0d2p-44L
 
-#define LN2O2    0x5.8b90bfbe8e7bcd6p-4L /* ln(2)/2 */
+#define LN2O2	0x5.8b90bfbe8e7bcd6p-4L /* ln(2)/2 */
+
+#define ZERO	0.L
 
-#define ZERO 	0.L
 /*
  * This set of coefficients is used in the polynomial approximation
  * for exp(r) where r is in [-ln(2)/2:ln(2)/2], and the r comes from
  * the argument reduction of x.
  */
-#define C0   0x1.000000000000000p0L
-#define C1   0x5.8b90bfbe8e7bcd6p-4L
-#define C2   0xf.5fdeffc162c7543p-8L
-#define C3   0x1.c6b08d704a0bf8bp-8L
-#define C4   0x2.76556df749cee54p-12L
-#define C5   0x2.bb0ffcf14ce6221p-16L
-#define C6   0x2.861225f0d8f0edfp-20L
-#define C7   0x1.ffcbfc588b0c687p-24L
-#define C8   0x1.62c0223a5c823fdp-28L
-#define C9   0xd.a929e9caf3e1ed2p-36L
-#define C10  0x7.933d4562e3b2cd7p-40L
-#define C11  0x3.d1958e6a3764b64p-44L
-#define C12  0x1.c3bd650fc1e343ap-48L
-#define C13  0xc.0b0c98b3649ff26p-56L
-#define C14  0x4.c525936609b02cfp-60L
-#define C15  0x1.c36e84400493e74p-64L
+#define C0	0x1.000000000000000p0L
+#define C1	0x5.8b90bfbe8e7bcd6p-4L
+#define C2	0xf.5fdeffc162c7543p-8L
+#define C3	0x1.c6b08d704a0bf8bp-8L
+#define C4	0x2.76556df749cee54p-12L
+#define C5	0x2.bb0ffcf14ce6221p-16L
+#define C6	0x2.861225f0d8f0edfp-20L
+#define C7	0x1.ffcbfc588b0c687p-24L
+#define C8	0x1.62c0223a5c823fdp-28L
+#define C9	0xd.a929e9caf3e1ed2p-36L
+#define C10	0x7.933d4562e3b2cd7p-40L
+#define C11	0x3.d1958e6a3764b64p-44L
+#define C12	0x1.c3bd650fc1e343ap-48L
+#define C13	0xc.0b0c98b3649ff26p-56L
+#define C14	0x4.c525936609b02cfp-60L
+#define C15	0x1.c36e84400493e74p-64L
 
 long double
 expl(long double x)
@@ -90,15 +91,10 @@
 	z.bits.sign = 0;
 
 	/* x is either 0 or a subnormal number. */
-	if (z.bits.exp == 0) {
-		if ((z.bits.manl | z.bits.manh) == 0)
-			return (1);
-		else
-			 return (1 + x);
-	}
+	if (z.bits.exp == 0)
+		return (1 + x);
 
 	if (XMIN <= x && x <= XMAX) {
-
 		/* Argument reduction. */
 		k = (int) (z.e * ILN2HI + z.e * ILN2LO);
 		r = z.e - k * LN2HI - k * LN2LO;
@@ -117,11 +113,13 @@
 		if (s) {
 			z.e = 1 / z.e;
 			z.bits.exp -= k;
-		} else
+		} else {
 			z.bits.exp += k;
+		}
 
 		return (z.e);
 	}
+
 	/*
 	 * If x = +Inf, then exp(x) = Inf.
 	 * If x = -Inf, then exp(x) = 0.
@@ -157,10 +155,11 @@
 	    (C6 + r * (C7 + r * (C8 + r * (C9 + r * (C10 + r * (C11 + r *
 	    (C12 + r * (C13 + r * (C14 + r * C15))))))))))))));
 	z.e = 1 / z.e;
+
 	/*
 	 * FIXME: There has to be a better way to handle gradual underflow
 	 * because the relative absolute error is fairly large for numerous
-	 * vlaues of x as exp(x) goes to zero.
+	 * values of x as exp(x) goes to zero.
 	 */
 	z.e = scalbnl(z.e, -k);
 


More information about the freebsd-standards mailing list