...
Run Format

Source file src/math/exp.go

Documentation: math

     1  // Copyright 2009 The Go Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package math
     6  
     7  // Exp returns e**x, the base-e exponential of x.
     8  //
     9  // Special cases are:
    10  //	Exp(+Inf) = +Inf
    11  //	Exp(NaN) = NaN
    12  // Very large values overflow to 0 or +Inf.
    13  // Very small values underflow to 1.
    14  func Exp(x float64) float64
    15  
    16  // The original C code, the long comment, and the constants
    17  // below are from FreeBSD's /usr/src/lib/msun/src/e_exp.c
    18  // and came with this notice. The go code is a simplified
    19  // version of the original C.
    20  //
    21  // ====================================================
    22  // Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
    23  //
    24  // Permission to use, copy, modify, and distribute this
    25  // software is freely granted, provided that this notice
    26  // is preserved.
    27  // ====================================================
    28  //
    29  //
    30  // exp(x)
    31  // Returns the exponential of x.
    32  //
    33  // Method
    34  //   1. Argument reduction:
    35  //      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
    36  //      Given x, find r and integer k such that
    37  //
    38  //               x = k*ln2 + r,  |r| <= 0.5*ln2.
    39  //
    40  //      Here r will be represented as r = hi-lo for better
    41  //      accuracy.
    42  //
    43  //   2. Approximation of exp(r) by a special rational function on
    44  //      the interval [0,0.34658]:
    45  //      Write
    46  //          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
    47  //      We use a special Remez algorithm on [0,0.34658] to generate
    48  //      a polynomial of degree 5 to approximate R. The maximum error
    49  //      of this polynomial approximation is bounded by 2**-59. In
    50  //      other words,
    51  //          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
    52  //      (where z=r*r, and the values of P1 to P5 are listed below)
    53  //      and
    54  //          |                  5          |     -59
    55  //          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
    56  //          |                             |
    57  //      The computation of exp(r) thus becomes
    58  //                             2*r
    59  //              exp(r) = 1 + -------
    60  //                            R - r
    61  //                                 r*R1(r)
    62  //                     = 1 + r + ----------- (for better accuracy)
    63  //                                2 - R1(r)
    64  //      where
    65  //                               2       4             10
    66  //              R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
    67  //
    68  //   3. Scale back to obtain exp(x):
    69  //      From step 1, we have
    70  //         exp(x) = 2**k * exp(r)
    71  //
    72  // Special cases:
    73  //      exp(INF) is INF, exp(NaN) is NaN;
    74  //      exp(-INF) is 0, and
    75  //      for finite argument, only exp(0)=1 is exact.
    76  //
    77  // Accuracy:
    78  //      according to an error analysis, the error is always less than
    79  //      1 ulp (unit in the last place).
    80  //
    81  // Misc. info.
    82  //      For IEEE double
    83  //          if x >  7.09782712893383973096e+02 then exp(x) overflow
    84  //          if x < -7.45133219101941108420e+02 then exp(x) underflow
    85  //
    86  // Constants:
    87  // The hexadecimal values are the intended ones for the following
    88  // constants. The decimal values may be used, provided that the
    89  // compiler will convert from decimal to binary accurately enough
    90  // to produce the hexadecimal values shown.
    91  
    92  func exp(x float64) float64 {
    93  	const (
    94  		Ln2Hi = 6.93147180369123816490e-01
    95  		Ln2Lo = 1.90821492927058770002e-10
    96  		Log2e = 1.44269504088896338700e+00
    97  
    98  		Overflow  = 7.09782712893383973096e+02
    99  		Underflow = -7.45133219101941108420e+02
   100  		NearZero  = 1.0 / (1 << 28) // 2**-28
   101  	)
   102  
   103  	// special cases
   104  	switch {
   105  	case IsNaN(x) || IsInf(x, 1):
   106  		return x
   107  	case IsInf(x, -1):
   108  		return 0
   109  	case x > Overflow:
   110  		return Inf(1)
   111  	case x < Underflow:
   112  		return 0
   113  	case -NearZero < x && x < NearZero:
   114  		return 1 + x
   115  	}
   116  
   117  	// reduce; computed as r = hi - lo for extra precision.
   118  	var k int
   119  	switch {
   120  	case x < 0:
   121  		k = int(Log2e*x - 0.5)
   122  	case x > 0:
   123  		k = int(Log2e*x + 0.5)
   124  	}
   125  	hi := x - float64(k)*Ln2Hi
   126  	lo := float64(k) * Ln2Lo
   127  
   128  	// compute
   129  	return expmulti(hi, lo, k)
   130  }
   131  
   132  // Exp2 returns 2**x, the base-2 exponential of x.
   133  //
   134  // Special cases are the same as Exp.
   135  func Exp2(x float64) float64
   136  
   137  func exp2(x float64) float64 {
   138  	const (
   139  		Ln2Hi = 6.93147180369123816490e-01
   140  		Ln2Lo = 1.90821492927058770002e-10
   141  
   142  		Overflow  = 1.0239999999999999e+03
   143  		Underflow = -1.0740e+03
   144  	)
   145  
   146  	// special cases
   147  	switch {
   148  	case IsNaN(x) || IsInf(x, 1):
   149  		return x
   150  	case IsInf(x, -1):
   151  		return 0
   152  	case x > Overflow:
   153  		return Inf(1)
   154  	case x < Underflow:
   155  		return 0
   156  	}
   157  
   158  	// argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2.
   159  	// computed as r = hi - lo for extra precision.
   160  	var k int
   161  	switch {
   162  	case x > 0:
   163  		k = int(x + 0.5)
   164  	case x < 0:
   165  		k = int(x - 0.5)
   166  	}
   167  	t := x - float64(k)
   168  	hi := t * Ln2Hi
   169  	lo := -t * Ln2Lo
   170  
   171  	// compute
   172  	return expmulti(hi, lo, k)
   173  }
   174  
   175  // exp1 returns e**r × 2**k where r = hi - lo and |r| ≤ ln(2)/2.
   176  func expmulti(hi, lo float64, k int) float64 {
   177  	const (
   178  		P1 = 1.66666666666666657415e-01  /* 0x3FC55555; 0x55555555 */
   179  		P2 = -2.77777777770155933842e-03 /* 0xBF66C16C; 0x16BEBD93 */
   180  		P3 = 6.61375632143793436117e-05  /* 0x3F11566A; 0xAF25DE2C */
   181  		P4 = -1.65339022054652515390e-06 /* 0xBEBBBD41; 0xC5D26BF1 */
   182  		P5 = 4.13813679705723846039e-08  /* 0x3E663769; 0x72BEA4D0 */
   183  	)
   184  
   185  	r := hi - lo
   186  	t := r * r
   187  	c := r - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))))
   188  	y := 1 - ((lo - (r*c)/(2-c)) - hi)
   189  	// TODO(rsc): make sure Ldexp can handle boundary k
   190  	return Ldexp(y, k)
   191  }
   192  

View as plain text