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