1*2f2c0062Sguenther /* $OpenBSD: e_expl.c,v 1.4 2016/09/12 19:47:02 guenther Exp $ */
249393c00Smartynas
349393c00Smartynas /*
449393c00Smartynas * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
549393c00Smartynas *
649393c00Smartynas * Permission to use, copy, modify, and distribute this software for any
749393c00Smartynas * purpose with or without fee is hereby granted, provided that the above
849393c00Smartynas * copyright notice and this permission notice appear in all copies.
949393c00Smartynas *
1049393c00Smartynas * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
1149393c00Smartynas * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
1249393c00Smartynas * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
1349393c00Smartynas * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
1449393c00Smartynas * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
1549393c00Smartynas * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
1649393c00Smartynas * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
1749393c00Smartynas */
1849393c00Smartynas
1949393c00Smartynas /* expl.c
2049393c00Smartynas *
2149393c00Smartynas * Exponential function, 128-bit long double precision
2249393c00Smartynas *
2349393c00Smartynas *
2449393c00Smartynas *
2549393c00Smartynas * SYNOPSIS:
2649393c00Smartynas *
2749393c00Smartynas * long double x, y, expl();
2849393c00Smartynas *
2949393c00Smartynas * y = expl( x );
3049393c00Smartynas *
3149393c00Smartynas *
3249393c00Smartynas *
3349393c00Smartynas * DESCRIPTION:
3449393c00Smartynas *
3549393c00Smartynas * Returns e (2.71828...) raised to the x power.
3649393c00Smartynas *
3749393c00Smartynas * Range reduction is accomplished by separating the argument
3849393c00Smartynas * into an integer k and fraction f such that
3949393c00Smartynas *
4049393c00Smartynas * x k f
4149393c00Smartynas * e = 2 e.
4249393c00Smartynas *
4349393c00Smartynas * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
4449393c00Smartynas * in the basic range [-0.5 ln 2, 0.5 ln 2].
4549393c00Smartynas *
4649393c00Smartynas *
4749393c00Smartynas * ACCURACY:
4849393c00Smartynas *
4949393c00Smartynas * Relative error:
5049393c00Smartynas * arithmetic domain # trials peak rms
5149393c00Smartynas * IEEE +-MAXLOG 100,000 2.6e-34 8.6e-35
5249393c00Smartynas *
5349393c00Smartynas *
5449393c00Smartynas * Error amplification in the exponential function can be
5549393c00Smartynas * a serious matter. The error propagation involves
5649393c00Smartynas * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
5749393c00Smartynas * which shows that a 1 lsb error in representing X produces
5849393c00Smartynas * a relative error of X times 1 lsb in the function.
5949393c00Smartynas * While the routine gives an accurate result for arguments
6049393c00Smartynas * that are exactly represented by a long double precision
6149393c00Smartynas * computer number, the result contains amplified roundoff
6249393c00Smartynas * error for large arguments not exactly represented.
6349393c00Smartynas *
6449393c00Smartynas *
6549393c00Smartynas * ERROR MESSAGES:
6649393c00Smartynas *
6749393c00Smartynas * message condition value returned
6849393c00Smartynas * exp underflow x < MINLOG 0.0
6949393c00Smartynas * exp overflow x > MAXLOG MAXNUM
7049393c00Smartynas *
7149393c00Smartynas */
7249393c00Smartynas
7349393c00Smartynas /* Exponential function */
7449393c00Smartynas
7549393c00Smartynas #include <float.h>
7649393c00Smartynas #include <math.h>
7749393c00Smartynas
7812129755Smartynas #include "math_private.h"
7912129755Smartynas
8049393c00Smartynas /* Pade' coefficients for exp(x) - 1
8149393c00Smartynas Theoretical peak relative error = 2.2e-37,
8249393c00Smartynas relative peak error spread = 9.2e-38
8349393c00Smartynas */
8449393c00Smartynas static long double P[5] = {
8549393c00Smartynas 3.279723985560247033712687707263393506266E-10L,
8649393c00Smartynas 6.141506007208645008909088812338454698548E-7L,
8749393c00Smartynas 2.708775201978218837374512615596512792224E-4L,
8849393c00Smartynas 3.508710990737834361215404761139478627390E-2L,
8949393c00Smartynas 9.999999999999999999999999999999999998502E-1L
9049393c00Smartynas };
9149393c00Smartynas static long double Q[6] = {
9249393c00Smartynas 2.980756652081995192255342779918052538681E-12L,
9349393c00Smartynas 1.771372078166251484503904874657985291164E-8L,
9449393c00Smartynas 1.504792651814944826817779302637284053660E-5L,
9549393c00Smartynas 3.611828913847589925056132680618007270344E-3L,
9649393c00Smartynas 2.368408864814233538909747618894558968880E-1L,
9749393c00Smartynas 2.000000000000000000000000000000000000150E0L
9849393c00Smartynas };
9949393c00Smartynas /* C1 + C2 = ln 2 */
1009fc13282Smartynas static const long double C1 = -6.93145751953125E-1L;
1019fc13282Smartynas static const long double C2 = -1.428606820309417232121458176568075500134E-6L;
10249393c00Smartynas
1039fc13282Smartynas static const long double LOG2EL = 1.442695040888963407359924681001892137426646L;
1049fc13282Smartynas static const long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L;
1059fc13282Smartynas static const long double MINLOGL = -1.143276959615573793352782661133116431383730e4L;
10649393c00Smartynas static const long double huge = 0x1p10000L;
10749393c00Smartynas #if 0 /* XXX Prevent gcc from erroneously constant folding this. */
10849393c00Smartynas static const long double twom10000 = 0x1p-10000L;
10949393c00Smartynas #else
11049393c00Smartynas static volatile long double twom10000 = 0x1p-10000L;
11149393c00Smartynas #endif
11249393c00Smartynas
11349393c00Smartynas long double
expl(long double x)11449393c00Smartynas expl(long double x)
11549393c00Smartynas {
11649393c00Smartynas long double px, xx;
11749393c00Smartynas int n;
11849393c00Smartynas
11949393c00Smartynas if( x > MAXLOGL)
12049393c00Smartynas return (huge*huge); /* overflow */
12149393c00Smartynas
12249393c00Smartynas if( x < MINLOGL )
12349393c00Smartynas return (twom10000*twom10000); /* underflow */
12449393c00Smartynas
12549393c00Smartynas /* Express e**x = e**g 2**n
12649393c00Smartynas * = e**g e**( n loge(2) )
12749393c00Smartynas * = e**( g + n loge(2) )
12849393c00Smartynas */
12949393c00Smartynas px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
13049393c00Smartynas n = px;
13149393c00Smartynas x += px * C1;
13249393c00Smartynas x += px * C2;
13349393c00Smartynas /* rational approximation for exponential
13449393c00Smartynas * of the fractional part:
13549393c00Smartynas * e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
13649393c00Smartynas */
13749393c00Smartynas xx = x * x;
13849393c00Smartynas px = x * __polevll( xx, P, 4 );
13949393c00Smartynas xx = __polevll( xx, Q, 5 );
14049393c00Smartynas x = px/( xx - px );
14149393c00Smartynas x = 1.0L + x + x;
14249393c00Smartynas
14349393c00Smartynas x = ldexpl( x, n );
14449393c00Smartynas return(x);
14549393c00Smartynas }
146*2f2c0062Sguenther DEF_STD(expl);
147