pow.js (5151B)
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 'use strict'; 20 21 // MODULES // 22 23 var isnan = require( './../../../../../base/assert/is-nan' ); 24 var PINF = require( '@stdlib/constants/float64/pinf' ); 25 26 27 // VARIABLES // 28 29 var ZERO = 0|0; // asm type annotation 30 var ONE = 1|0; // asm type annotation 31 32 33 // MAIN // 34 35 /** 36 * Evaluates the exponential function. 37 * 38 * ## Method 39 * 40 * - The exponential function is given by 41 * 42 * ```tex 43 * z = x^y 44 * ``` 45 * 46 * where \\(x\\) is the base and \\(y\\) the exponent. 47 * 48 * - First observe that a naive approach for exponentiation 49 * 50 * ```tex 51 * 5^5 = 5 \cdot 5 \cdot 5 \cdot 5 \cdot 5 52 * ``` 53 * 54 * requires \\(y-1\\) multiplications. 55 * 56 * - We can reduce the number of multiplications by first computing \\(x2 = x \cdot x\\). 57 * 58 * ```tex 59 * 5^5 = x2 \cdot x2 \cdot x 60 * ``` 61 * 62 * thus requiring only three multiplications. 63 * 64 * - This observation may be generalized, such that, for a positive exponent \\(y\\), 65 * 66 * ```tex 67 * x^y = \begin{cases} 68 * x (x^2)^{\frac{y-1}{2}}, & \text{if $y$ is odd} \\ 69 * (x^2)^{\frac{y}{2}}, & \text{if $y$ is even} 70 * \end{cases} 71 * ``` 72 * 73 * - Note that the above generalization only involves powers of two. For example, in our working example, the powers are \\(1\\) and \\(4\\). To determine these powers, we observe that integer values, when stored in binary format, are simply sums of powers of two. For example, the integer \\(5\\) has the bit sequence 74 * 75 * ```binarystring 76 * 00000000000000000000000000000101 77 * ``` 78 * 79 * where \\(101\\) translates to 80 * 81 * ```tex 82 * 2^2 + 2^0 = 4 + 1 = 5 83 * ``` 84 * 85 * Thus, rather conveniently, the powers of two needed for exponentiation are easily derived from the binary representation of the integer exponent. 86 * 87 * - The previous observation lends itself readily to an iterative exponentiation algorithm, referred to as **right-to-left binary exponentiation**. The algorithm is as follows: 88 * 89 * ```text 90 * 1. Examine the least significant bit to determine if we have a power of 2. 91 * 2. If yes, compute an intermediate result. 92 * 3. Square the base. 93 * 4. Shift off the least significant bit (LSB). 94 * 5. If the exponent is greater than 0, repeat steps 1-4. 95 * 6. Return the intermediate result. 96 * ``` 97 * 98 * For example, consider \\(5^5 = 3125\\). 99 * 100 * ```text 101 * Initialization: r = 1 102 * Iteration 1: y = 101 => r = 1*5, x = 5*5 = 25 103 * Iteration 2: y = 10 => x = 25*25 = 625 104 * Iteration 3: y = 1 => r = 5*625 = 3125, x = 625*625 105 * Return: r 106 * ``` 107 * 108 * 109 * ## Notes 110 * 111 * - The above algorithm involves \\(\lfloor \log_2(y) \rfloor\\) square operations and at most \\(\lfloor \log_2(y) \rfloor\\) multiplications. 112 * 113 * - The above algorithm may not return precise results due to an accumulation of error. For example, 114 * 115 * ```javascript 116 * var y = pow( 10.0, 308 ); 117 * // returns 1.0000000000000006e+308 118 * // expected 1.0e+308 119 * ``` 120 * 121 * If we compare the bit sequence of the returned value 122 * 123 * ```binarystring 124 * 0111111111100001110011001111001110000101111010111100100010100011 125 * ``` 126 * 127 * with the expected value 128 * 129 * ```binarystring 130 * 0111111111100001110011001111001110000101111010111100100010100000 131 * ``` 132 * 133 * we observe that the returned value differs in its last two bits. 134 * 135 * 136 * @param {number} x - base 137 * @param {integer32} y - exponent 138 * @returns {number} function value 139 * 140 * @example 141 * var v = pow( 2.0, 3 ); 142 * // returns 8.0 143 * 144 * @example 145 * var v = pow( 3.14, 0 ); 146 * // returns 1.0 147 * 148 * @example 149 * var v = pow( 2.0, -2 ); 150 * // returns 0.25 151 * 152 * @example 153 * var v = pow( 0.0, 0 ); 154 * // returns 1.0 155 * 156 * @example 157 * var v = pow( -3.14, 1 ); 158 * // returns -3.14 159 * 160 * @example 161 * var v = pow( NaN, 0 ); 162 * // returns NaN 163 */ 164 function pow( x, y ) { 165 var v; 166 167 if ( isnan( x ) ) { 168 return NaN; 169 } 170 // If the exponent is negative, use the reciprocal... 171 if ( y < ZERO ) { 172 y = -y; 173 if ( x === 0.0 ) { 174 x = 1.0 / x; // +-infinity 175 if ( ( y & ONE ) === ONE ) { 176 // Exponent is odd, so `x` keeps its sign: 177 return x; 178 } 179 // Exponent is even, so result is always positive: 180 return PINF; 181 } 182 x = 1.0 / x; 183 } 184 // If the exponent is zero, the result is always unity... 185 else if ( y === ZERO ) { 186 return 1.0; 187 } 188 v = 1; 189 while ( y !== ZERO ) { 190 // Check the least significant bit (LSB) to determine if "on" (if so, we have a power of 2)... 191 if ( ( y & ONE ) === ONE ) { 192 v *= x; 193 } 194 x *= x; // possible overflow 195 y >>= ONE; 196 } 197 return v; 198 } 199 200 201 // EXPORTS // 202 203 module.exports = pow;