Implementation of expl()

Sun Dec 9 13:39:17 PST 2007

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.

I don't remember if bde sent me comments on this code.  I sure
he has plenty. :)  

steve

/*-
 * Copyright (c) 2007 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * Use argument reduction to compute exp(x).  The reduction writes
 * x = k * ln(2) + r with r in the range [0, ln(2)].  This then
 * gives exp(x) = 2**k * exp(r), and exp(r) is evaluated via a nearly
 * minimax polynomial approximation such that r is mapped into [-1:1].
 */
#include "math.h"
#include "math_private.h"
#include "fpmath.h"

/* ln(LDBL_MAX) =  11356.523406294144 */
#define XMAX	0x2.c5c85fdf473de6ap12L

/* ln(LDBL_MIN) = -11355.137111933024 */
#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 TOLN2    0x2.e2a8eca5705fc2fp0L  /* 2/ln(2) */

#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 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

long double
expl(long double x)
{
	union IEEEl2bits z;
	int k, s;
	long double r;

	z.e = x;
	s = z.bits.sign;
	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 (XMIN <= x && x <= XMAX) {

		/* Argument reduction. */
		k = (int) (z.e * ILN2HI + z.e * ILN2LO);
		r = z.e - k * LN2HI - k * LN2LO;
		if (r > LN2O2) {
			r -= LN2;
			k++;
		}

		/* Compute exp(r) via the polynomial approximation. */
		r = r * TOLN2;
		z.e = C0 + r * (C1 + r * (C2 + r * (C3 + r * (C4 + r *
		    (C5 + r * (C6 + r * (C7 + r * (C8 + r * (C9 + r *
		    (C10 + r * (C11 + r * (C12 + r * (C13 + r *
		    (C14 + r * C15))))))))))))));

		if (s) {
			z.e = 1 / z.e;
			z.bits.exp -= k;
		} else
			z.bits.exp += k;

		return (z.e);
	}
	/*
	 * If x = +Inf, then exp(x) = Inf.
	 * If x = -Inf, then exp(x) = 0.
	 * If x = NaN, then exp(x) = NaN.
	 */
	if (z.bits.exp == 32767) {
		mask_nbit_l(z);
		if (!(z.bits.manh | z.bits.manl))
			return (s ? ZERO : 1.L / ZERO);
		return ((x - x) / (x - x));
	}

	/* If x > 0 then, overflow to +Inf. */
	if (!s)
		return (1.L / ZERO);

	/* For x < 0, check if gradual underflow is needed. */
	if (z.e > GRAD)
		return (ZERO);

	/* Argument reduction. */
	k = (int) (z.e * ILN2);
	r = z.e - k * LN2HI - k * LN2LO;
	if (r > LN2O2) {
		r -= LN2;
		k++;
	}

	/* Compute exp(r) via the polynomial approximation. */
	r *= TOLN2;

	z.e = C0 + r * (C1 + r * (C2 + r * (C3 + r * (C4 + r * (C5 + r *
	    (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.
	 */
	z.e = scalbnl(z.e, -k);

	return (z.e);
}