xref: /openbsd-src/lib/libm/src/ld128/e_expl.c (revision 2f2c00629eff6a304ebffb255fc56f4fa7a1833b)
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