compiler/stdc/math/ld128/e_expl.c

   1 /*      $OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:18 martynas Exp $     */
   2
   3 /*
   4  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
   5  *
   6  * Permission to use, copy, modify, and distribute this software for any
   7  * purpose with or without fee is hereby granted, provided that the above
   8  * copyright notice and this permission notice appear in all copies.
   9  *
  10  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  11  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  12  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  13  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  14  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  15  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  16  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  17  */
  18
  19 /*                                                      expl.c
  20  *
  21  *      Exponential function, 128-bit long double precision
  22  *
  23  *
  24  *
  25  * SYNOPSIS:
  26  *
  27  * long double x, y, expl();
  28  *
  29  * y = expl( x );
  30  *
  31  *
  32  *
  33  * DESCRIPTION:
  34  *
  35  * Returns e (2.71828...) raised to the x power.
  36  *
  37  * Range reduction is accomplished by separating the argument
  38  * into an integer k and fraction f such that
  39  *
  40  *     x    k  f
  41  *    e  = 2  e.
  42  *
  43  * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
  44  * in the basic range [-0.5 ln 2, 0.5 ln 2].
  45  *
  46  *
  47  * ACCURACY:
  48  *
  49  *                      Relative error:
  50  * arithmetic   domain     # trials      peak         rms
  51  *    IEEE      +-MAXLOG    100,000     2.6e-34     8.6e-35
  52  *
  53  *
  54  * Error amplification in the exponential function can be
  55  * a serious matter.  The error propagation involves
  56  * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
  57  * which shows that a 1 lsb error in representing X produces
  58  * a relative error of X times 1 lsb in the function.
  59  * While the routine gives an accurate result for arguments
  60  * that are exactly represented by a long double precision
  61  * computer number, the result contains amplified roundoff
  62  * error for large arguments not exactly represented.
  63  *
  64  *
  65  * ERROR MESSAGES:
  66  *
  67  *   message         condition      value returned
  68  * exp underflow    x < MINLOG         0.0
  69  * exp overflow     x > MAXLOG         MAXNUM
  70  *
  71  */
  72
  73 /*      Exponential function    */
  74
  75 #include <float.h>
  76 #include "math.h"
  77
  78 #include "math_private.h"
  79
  80 /* Pade' coefficients for exp(x) - 1
  81    Theoretical peak relative error = 2.2e-37,
  82    relative peak error spread = 9.2e-38
  83  */
  84 static long double P[5] = {
  85  3.279723985560247033712687707263393506266E-10L,
  86  6.141506007208645008909088812338454698548E-7L,
  87  2.708775201978218837374512615596512792224E-4L,
  88  3.508710990737834361215404761139478627390E-2L,
  89  9.999999999999999999999999999999999998502E-1L
  90 };
  91 static long double Q[6] = {
  92  2.980756652081995192255342779918052538681E-12L,
  93  1.771372078166251484503904874657985291164E-8L,
  94  1.504792651814944826817779302637284053660E-5L,
  95  3.611828913847589925056132680618007270344E-3L,
  96  2.368408864814233538909747618894558968880E-1L,
  97  2.000000000000000000000000000000000000150E0L
  98 };
  99 /* C1 + C2 = ln 2 */
 100 static const long double C1 = -6.93145751953125E-1L;
 101 static const long double C2 = -1.428606820309417232121458176568075500134E-6L;
 102
 103 static const long double LOG2EL = 1.442695040888963407359924681001892137426646L;
 104 static const long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L;
 105 static const long double MINLOGL = -1.143276959615573793352782661133116431383730e4L;
 106 static const long double huge = 0x1p10000L;
 107 #if 0 /* XXX Prevent gcc from erroneously constant folding this. */
 108 static const long double twom10000 = 0x1p-10000L;
 109 #else
 110 static const volatile long double twom10000 __attribute__ ((__section__(".rodata,\"a\" " SECTIONCOMMENT))) = 0x1p-10000L;
 111 #endif
 112
 113 long double
 114 expl(long double x)
 115 {
 116 long double px, xx;
 117 int n;
 118
 119 if( x > MAXLOGL)
 120         return (huge*huge);             /* overflow */
 121
 122 if( x < MINLOGL )
 123         return (twom10000*twom10000);   /* underflow */
 124
 125 /* Express e**x = e**g 2**n
 126  *   = e**g e**( n loge(2) )
 127  *   = e**( g + n loge(2) )
 128  */
 129 px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
 130 n = px;
 131 x += px * C1;
 132 x += px * C2;
 133 /* rational approximation for exponential
 134  * of the fractional part:
 135  * e**x =  1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
 136  */
 137 xx = x * x;
 138 px = x * __polevll( xx, P, 4 );
 139 xx = __polevll( xx, Q, 5 );
 140 x =  px/( xx - px );
 141 x = 1.0L + x + x;
 142
 143 x = ldexpl( x, n );
 144 return(x);
 145 }