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