exp.js (5298B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2018 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 * 18 * 19 * ## Notice 20 * 21 * The following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/e_exp.c}. The implementation follows the original, but has been modified for JavaScript. 22 * 23 * ```text 24 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. 25 * 26 * Developed at SunPro, a Sun Microsystems, Inc. business. 27 * Permission to use, copy, modify, and distribute this 28 * software is freely granted, provided that this notice 29 * is preserved. 30 * ``` 31 */ 32 33 'use strict'; 34 35 // MODULES // 36 37 var isnan = require( './../../../../base/assert/is-nan' ); 38 var trunc = require( './../../../../base/special/trunc' ); 39 var NINF = require( '@stdlib/constants/float64/ninf' ); 40 var PINF = require( '@stdlib/constants/float64/pinf' ); 41 var expmulti = require( './expmulti.js' ); 42 43 44 // VARIABLES // 45 46 var LN2_HI = 6.93147180369123816490e-01; 47 var LN2_LO = 1.90821492927058770002e-10; 48 var LOG2_E = 1.44269504088896338700e+00; 49 var OVERFLOW = 7.09782712893383973096e+02; 50 var UNDERFLOW = -7.45133219101941108420e+02; 51 var NEARZERO = 1.0 / (1 << 28); // 2^-28; 52 var NEG_NEARZERO = -NEARZERO; 53 54 55 // MAIN // 56 57 /** 58 * Evaluates the natural exponential function. 59 * 60 * ## Method 61 * 62 * 1. We reduce \\( x \\) to an \\( r \\) so that \\( |r| \leq 0.5 \cdot \ln(2) \approx 0.34658 \\). Given \\( x \\), we find an \\( r \\) and integer \\( k \\) such that 63 * 64 * ```tex 65 * \begin{align*} 66 * x &= k \cdot \ln(2) + r \\ 67 * |r| &\leq 0.5 \cdot \ln(2) 68 * \end{align*} 69 * ``` 70 * 71 * <!-- <note> --> 72 * 73 * \\( r \\) can be represented as \\( r = \mathrm{hi} - \mathrm{lo} \\) for better accuracy. 74 * 75 * <!-- </note> --> 76 * 77 * 2. We approximate of \\( e^{r} \\) by a special rational function on the interval \\(\[0,0.34658]\\): 78 * 79 * ```tex 80 * \begin{align*} 81 * R\left(r^2\right) &= r \cdot \frac{ e^{r}+1 }{ e^{r}-1 } \\ 82 * &= 2 + \frac{r^2}{6} - \frac{r^4}{360} + \ldots 83 * \end{align*} 84 * ``` 85 * 86 * We use a special Remes algorithm on \\(\[0,0.34658]\\) to generate a polynomial of degree \\(5\\) to approximate \\(R\\). The maximum error of this polynomial approximation is bounded by \\(2^{-59}\\). In other words, 87 * 88 * ```tex 89 * R(z) \sim 2 + P_1 z + P_2 z^2 + P_3 z^3 + P_4 z^4 + P_5 z^5 90 * ``` 91 * 92 * where \\( z = r^2 \\) and 93 * 94 * ```tex 95 * \left| 2 + P_1 z + \ldots + P_5 z^5 - R(z) \right| \leq 2^{-59} 96 * ``` 97 * 98 * <!-- <note> --> 99 * 100 * The values of \\( P_1 \\) to \\( P_5 \\) are listed in the source code. 101 * 102 * <!-- </note> --> 103 * 104 * The computation of \\( e^{r} \\) thus becomes 105 * 106 * ```tex 107 * \begin{align*} 108 * e^{r} &= 1 + \frac{2r}{R-r} \\ 109 * &= 1 + r + \frac{r \cdot R_1(r)}{2 - R_1(r)}\ \text{for better accuracy} 110 * \end{align*} 111 * ``` 112 * 113 * where 114 * 115 * ```tex 116 * R_1(r) = r - P_1\ r^2 + P_2\ r^4 + \ldots + P_5\ r^{10} 117 * ``` 118 * 119 * 3. We scale back to obtain \\( e^{x} \\). From step 1, we have 120 * 121 * ```tex 122 * e^{x} = 2^k e^{r} 123 * ``` 124 * 125 * 126 * ## Special Cases 127 * 128 * ```tex 129 * \begin{align*} 130 * e^\infty &= \infty \\ 131 * e^{-\infty} &= 0 \\ 132 * e^{\mathrm{NaN}} &= \mathrm{NaN} \\ 133 * e^0 &= 1\ \mathrm{is\ exact\ for\ finite\ argument\ only} 134 * \end{align*} 135 * ``` 136 * 137 * ## Notes 138 * 139 * - According to an error analysis, the error is always less than \\(1\\) ulp (unit in the last place). 140 * 141 * - For an IEEE double, 142 * 143 * - if \\(x > 7.09782712893383973096\mbox{e+}02\\), then \\(e^{x}\\) overflows 144 * - if \\(x < -7.45133219101941108420\mbox{e+}02\\), then \\(e^{x}\\) underflows 145 * 146 * - The hexadecimal values included in the source code are the intended ones for the used constants. Decimal values may be used, provided that the compiler will convert from decimal to binary accurately enough to produce the intended hexadecimal values. 147 * 148 * 149 * @param {number} x - input value 150 * @returns {number} function value 151 * 152 * @example 153 * var v = exp( 4.0 ); 154 * // returns ~54.5982 155 * 156 * @example 157 * var v = exp( -9.0 ); 158 * // returns ~1.234e-4 159 * 160 * @example 161 * var v = exp( 0.0 ); 162 * // returns 1.0 163 * 164 * @example 165 * var v = exp( NaN ); 166 * // returns NaN 167 */ 168 function exp( x ) { 169 var hi; 170 var lo; 171 var k; 172 173 if ( isnan( x ) || x === PINF ) { 174 return x; 175 } 176 if ( x === NINF ) { 177 return 0.0; 178 } 179 if ( x > OVERFLOW ) { 180 return PINF; 181 } 182 if ( x < UNDERFLOW ) { 183 return 0.0; 184 } 185 if ( 186 x > NEG_NEARZERO && 187 x < NEARZERO 188 ) { 189 return 1.0 + x; 190 } 191 // Reduce and compute `r = hi - lo` for extra precision. 192 if ( x < 0.0 ) { 193 k = trunc( (LOG2_E*x) - 0.5 ); 194 } else { 195 k = trunc( (LOG2_E*x) + 0.5 ); 196 } 197 hi = x - (k*LN2_HI); 198 lo = k * LN2_LO; 199 200 return expmulti( hi, lo, k ); 201 } 202 203 204 // EXPORTS // 205 206 module.exports = exp;