log2ax.js (5904B)
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 and license were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/9.3.0/lib/msun/src/s_pow.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 getHighWord = require( '@stdlib/number/float64/base/get-high-word' ); 38 var setLowWord = require( '@stdlib/number/float64/base/set-low-word' ); 39 var setHighWord = require( '@stdlib/number/float64/base/set-high-word' ); 40 var BIAS = require( '@stdlib/constants/float64/exponent-bias' ); 41 var polyvalL = require( './polyval_l.js' ); 42 43 44 // VARIABLES // 45 46 // 0x000fffff = 1048575 => 0 00000000000 11111111111111111111 47 var HIGH_SIGNIFICAND_MASK = 0x000fffff|0; // asm type annotation 48 49 // 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022 50 var HIGH_MIN_NORMAL_EXP = 0x00100000|0; // asm type annotation 51 52 // 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1 53 var HIGH_BIASED_EXP_0 = 0x3ff00000|0; // asm type annotation 54 55 // 0x20000000 = 536870912 => 0 01000000000 00000000000000000000 => biased exponent: 512 = -511+1023 56 var HIGH_BIASED_EXP_NEG_512 = 0x20000000|0; // asm type annotation 57 58 // 0x00080000 = 524288 => 0 00000000000 10000000000000000000 59 var HIGH_SIGNIFICAND_HALF = 0x00080000|0; // asm type annotation 60 61 // TODO: consider making an external constant 62 var HIGH_NUM_SIGNIFICAND_BITS = 20|0; // asm type annotation 63 64 var TWO53 = 9007199254740992.0; // 0x43400000, 0x00000000 65 66 // 2/(3*LN2) 67 var CP = 9.61796693925975554329e-01; // 0x3FEEC709, 0xDC3A03FD 68 69 // (float)CP 70 var CP_HI = 9.61796700954437255859e-01; // 0x3FEEC709, 0xE0000000 71 72 // Low: CP_HI 73 var CP_LO = -7.02846165095275826516e-09; // 0xBE3E2FE0, 0x145B01F5 74 75 var BP = [ 76 1.0, 77 1.5 78 ]; 79 var DP_HI = [ 80 0.0, 81 5.84962487220764160156e-01 // 0x3FE2B803, 0x40000000 82 ]; 83 var DP_LO = [ 84 0.0, 85 1.35003920212974897128e-08 // 0x3E4CFDEB, 0x43CFD006 86 ]; 87 88 89 // MAIN // 90 91 /** 92 * Computes \\(\operatorname{log2}(ax)\\). 93 * 94 * @private 95 * @param {Array} out - output array 96 * @param {number} ax - absolute value of `x` 97 * @param {number} ahx - high word of `ax` 98 * @returns {Array} output array containing a tuple comprised of high and low parts 99 * 100 * @example 101 * var t = log2ax( [ 0.0, 0.0 ], 9.0, 1075970048 ); // => [ t1, t2 ] 102 * // returns [ 3.169923782348633, 0.0000012190936795504075 ] 103 */ 104 function log2ax( out, ax, ahx ) { 105 var tmp; 106 var ss; // `hs + ls` 107 var s2; // `ss` squared 108 var hs; 109 var ls; 110 var ht; 111 var lt; 112 var bp; // `BP` constant 113 var dp; // `DP` constant 114 var hp; 115 var lp; 116 var hz; 117 var lz; 118 var t1; 119 var t2; 120 var t; 121 var r; 122 var u; 123 var v; 124 var n; 125 var j; 126 var k; 127 128 n = 0|0; // asm type annotation 129 130 // Check if `x` is subnormal... 131 if ( ahx < HIGH_MIN_NORMAL_EXP ) { 132 ax *= TWO53; 133 n -= 53|0; // asm type annotation 134 ahx = getHighWord( ax ); 135 } 136 // Extract the unbiased exponent of `x`: 137 n += ((ahx >> HIGH_NUM_SIGNIFICAND_BITS) - BIAS)|0; // asm type annotation 138 139 // Isolate the significand bits of `x`: 140 j = (ahx & HIGH_SIGNIFICAND_MASK)|0; // asm type annotation 141 142 // Normalize `ahx` by setting the (biased) exponent to `1023`: 143 ahx = (j | HIGH_BIASED_EXP_0)|0; // asm type annotation 144 145 // Determine the interval of `|x|` by comparing significand bits... 146 147 // |x| < sqrt(3/2) 148 if ( j <= 0x3988E ) { // 0 00000000000 00111001100010001110 149 k = 0; 150 } 151 // |x| < sqrt(3) 152 else if ( j < 0xBB67A ) { // 0 00000000000 10111011011001111010 153 k = 1; 154 } 155 // |x| >= sqrt(3) 156 else { 157 k = 0; 158 n += 1|0; // asm type annotation 159 ahx -= HIGH_MIN_NORMAL_EXP; 160 } 161 // Load the normalized high word into `|x|`: 162 ax = setHighWord( ax, ahx ); 163 164 // Compute `ss = hs + ls = (x-1)/(x+1)` or `(x-1.5)/(x+1.5)`: 165 bp = BP[ k ]; // BP[0] = 1.0, BP[1] = 1.5 166 u = ax - bp; // (x-1) || (x-1.5) 167 v = 1.0 / (ax + bp); // 1/(x+1) || 1/(x+1.5) 168 ss = u * v; 169 hs = setLowWord( ss, 0 ); // set all low word (less significant significand) bits to 0s 170 171 // Compute `ht = ax + bp` (via manipulation, i.e., bit flipping, of the high word): 172 tmp = ((ahx>>1) | HIGH_BIASED_EXP_NEG_512) + HIGH_SIGNIFICAND_HALF; 173 tmp += (k << 18); // `(k<<18)` can be considered the word equivalent of `1.0` or `1.5` 174 ht = setHighWord( 0.0, tmp ); 175 lt = ax - (ht - bp); 176 ls = v * ( ( u - (hs*ht) ) - ( hs*lt ) ); 177 178 // Compute `log(ax)`... 179 180 s2 = ss * ss; 181 r = s2 * s2 * polyvalL( s2 ); 182 r += ls * (hs + ss); 183 s2 = hs * hs; 184 ht = 3.0 + s2 + r; 185 ht = setLowWord( ht, 0 ); 186 lt = r - ((ht-3.0) - s2); 187 188 // u+v = ss*(1+...): 189 u = hs * ht; 190 v = ( ls*ht ) + ( lt*ss ); 191 192 // 2/(3LN2) * (ss+...): 193 hp = u + v; 194 hp = setLowWord( hp, 0 ); 195 lp = v - (hp - u); 196 hz = CP_HI * hp; // CP_HI+CP_LO = 2/(3*LN2) 197 lz = ( CP_LO*hp ) + ( lp*CP ) + DP_LO[ k ]; 198 199 // log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz 200 dp = DP_HI[ k ]; 201 t = n; 202 t1 = ((hz+lz) + dp) + t; // log2(ax) 203 t1 = setLowWord( t1, 0 ); 204 t2 = lz - (((t1-t) - dp) - hz); 205 206 out[ 0 ] = t1; 207 out[ 1 ] = t2; 208 return out; 209 } 210 211 212 // EXPORTS // 213 214 module.exports = log2ax;