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