|
| 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/e_powf.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 fromWordf = require( '@stdlib/number/float32/base/from-word' ); |
| 38 | +var toWordf = require( '@stdlib/number/float32/base/to-word' ); |
| 39 | +var float64ToFloat32 = require( '@stdlib/number/float64/base/to-float32' ); |
| 40 | +var BIAS = require( '@stdlib/constants/float32/exponent-bias' ); |
| 41 | +var polyvalL = require( './polyval_l.js' ); |
| 42 | + |
| 43 | + |
| 44 | +// VARIABLES // |
| 45 | + |
| 46 | +// 0x007fffff = 8388607 => 0 00000000000 11111111111111111111 |
| 47 | +var HIGH_EXP_2 = 0x007fffff; // asm type annotation |
| 48 | + |
| 49 | +// 0x00080000 = 524288 => 0 00000000000 10000000000000000000 |
| 50 | +var HIGH_MIN_NORMAL_EXP = 0x00800000; // asm type annotation |
| 51 | + |
| 52 | +// 0x3f800000 = 1065353216 => 0 01111111000 00000000000000000000 |
| 53 | +var HIGH_BIASED_EXP_0 = 0x3f800000; // asm type annotation |
| 54 | + |
| 55 | +// 0xfffff000 = 4294963200 => 1 11111111111 11111111111100000000 |
| 56 | +var HIGH_BIASED_EXP_NEG_64 = 0xfffff000; // asm type annotation |
| 57 | + |
| 58 | +// 0x20000000 = 536870912 => 0 10000000000 00000000000000000000 |
| 59 | +var HIGH_SIGNIFICAND_HALF = 0x20000000; // asm type annotation |
| 60 | + |
| 61 | +// 0x00400000 = 4194304 => 0 00000000000 10000000000000000000 |
| 62 | +var HIGH_MIN_NORMAL_EXP_HALF = 0x00400000; // asm type annotation |
| 63 | + |
| 64 | +var HIGH_NUM_SIGNIFICAND_BITS = 23; // asm type annotation |
| 65 | + |
| 66 | +var TWO24 = 16777216.0; // 0x4b800000 |
| 67 | + |
| 68 | +// 2/(3*LN2) |
| 69 | +var CP = 9.6179670095e-01; // 0x3f76384f =2/(3ln2) |
| 70 | + |
| 71 | +// (float)CP |
| 72 | +var CP_HI = 9.6191406250e-01; // 0x3f764000 =12b cp |
| 73 | + |
| 74 | +// Low: CP_HI |
| 75 | +var CP_LO = -1.1736857402e-04; // 0xb8f623c6 =tail of CP_HI |
| 76 | + |
| 77 | +var BP = [ |
| 78 | + 1.0, |
| 79 | + 1.5 |
| 80 | +]; |
| 81 | +var DP_HI = [ |
| 82 | + 0.0, |
| 83 | + 5.84960938e-01 // 0x3f15c000 |
| 84 | +]; |
| 85 | +var DP_LO = [ |
| 86 | + 0.0, |
| 87 | + 1.56322085e-06 // 0x35d1cfdc |
| 88 | +]; |
| 89 | + |
| 90 | + |
| 91 | +// MAIN // |
| 92 | + |
| 93 | +/** |
| 94 | +* Computes \\(\operatorname{log2}(ax)\\). |
| 95 | +* |
| 96 | +* @private |
| 97 | +* @param {Array} out - output array |
| 98 | +* @param {number} ax - absolute value of `x` |
| 99 | +* @param {number} ahx - high word of `ax` |
| 100 | +* @returns {Array} output array containing a tuple comprised of high and low parts |
| 101 | +* |
| 102 | +* @example |
| 103 | +* var t = log2ax( [ 0.0, 0.0 ], 9.0, 1075970048 ); // => [ t1, t2 ] |
| 104 | +* // returns [ 3.169923782348633, 0.0000012190936795504075 ] |
| 105 | +*/ |
| 106 | +function log2ax( out, ax, ahx ) { |
| 107 | + var ahs; |
| 108 | + var ss; // `hs + ls` |
| 109 | + var s2; // `ss` squared |
| 110 | + var hs; |
| 111 | + var ls; |
| 112 | + var ht; |
| 113 | + var lt; |
| 114 | + var bp; // `BP` constant |
| 115 | + var dp; // `DP` constant |
| 116 | + var hp; |
| 117 | + var lp; |
| 118 | + var hz; |
| 119 | + var lz; |
| 120 | + var t1; |
| 121 | + var t2; |
| 122 | + var t; |
| 123 | + var r; |
| 124 | + var u; |
| 125 | + var v; |
| 126 | + var n; |
| 127 | + var j; |
| 128 | + var k; |
| 129 | + |
| 130 | + n = 0; // asm type annotation |
| 131 | + |
| 132 | + // Check if `x` is subnormal... |
| 133 | + if ( ahx < HIGH_MIN_NORMAL_EXP ) { |
| 134 | + ax *= TWO24; |
| 135 | + n -= 24; // asm type annotation |
| 136 | + ahx = fromWordf( ax ); |
| 137 | + } |
| 138 | + // Extract the unbiased exponent of `x`: |
| 139 | + n += ((ahx >> HIGH_NUM_SIGNIFICAND_BITS) - BIAS); // asm type annotation |
| 140 | + |
| 141 | + // Isolate the significand bits of `x`: |
| 142 | + j = (ahx & HIGH_EXP_2); // asm type annotation |
| 143 | + |
| 144 | + // Normalize `ahx` by setting the (biased) exponent to `127`: |
| 145 | + ahx = (j | HIGH_BIASED_EXP_0); // asm type annotation |
| 146 | + |
| 147 | + // Determine the interval of `|x|` by comparing significand bits... |
| 148 | + |
| 149 | + // |x| < sqrt(3/2) |
| 150 | + if ( j <= 0x1cc471 ) { // 0 00000000001 11001100100010001110 |
| 151 | + k = 0; |
| 152 | + } |
| 153 | + // |x| < sqrt(3) |
| 154 | + else if ( j < 0x5db3d7 ) { // 0 00000000010 11101101100111110111 |
| 155 | + k = 1; |
| 156 | + } |
| 157 | + // |x| >= sqrtf(3) |
| 158 | + else { |
| 159 | + k = 0; |
| 160 | + n += 1; // asm type annotation |
| 161 | + ahx -= HIGH_MIN_NORMAL_EXP; |
| 162 | + } |
| 163 | + // Load the normalized float word into `|x|`: |
| 164 | + ahx = float64ToFloat32( ahx ); |
| 165 | + ax = toWordf( ahx ); |
| 166 | + |
| 167 | + // Compute `ss = hs + ls = (x-1)/(x+1)` or `(x-1.5)/(x+1.5)`: |
| 168 | + bp = float64ToFloat32( BP[ k ] ); // BP[0] = 1.0, BP[1] = 1.5 |
| 169 | + u = float64ToFloat32( ax - bp ); // (x-1) || (x-1.5) |
| 170 | + v = float64ToFloat32( 1.0 / float64ToFloat32(ax + bp) ); // 1/(x+1) || 1/(x+1.5) |
| 171 | + ss = float64ToFloat32( u * v ); |
| 172 | + hs = float64ToFloat32( ss ); |
| 173 | + ahs = fromWordf( hs ); |
| 174 | + ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ) |
| 175 | + hs = toWordf( ahs ); |
| 176 | + |
| 177 | + // Compute `ht = ax + bp` High |
| 178 | + ahs = ((ahx>>1) & HIGH_BIASED_EXP_NEG_64) | HIGH_SIGNIFICAND_HALF; |
| 179 | + ahs += HIGH_MIN_NORMAL_EXP_HALF+(k << 21); // `(k<<21)` can be considered the word equivalent of `1.0` or `1.5` |
| 180 | + ahs = float64ToFloat32( ahs ); |
| 181 | + ht = toWordf( ahs ); |
| 182 | + lt = float64ToFloat32( ax - float64ToFloat32(ht - bp) ); |
| 183 | + ls = float64ToFloat32( v * ( float64ToFloat32( u - float64ToFloat32(hs*ht) ) - float64ToFloat32( hs*lt ) ) ); // eslint-disable-line max-len |
| 184 | + |
| 185 | + // Compute `log(ax)`... |
| 186 | + |
| 187 | + s2 = float64ToFloat32( ss * ss); |
| 188 | + r = float64ToFloat32( s2 * s2 * float64ToFloat32( polyvalL( s2 ))); |
| 189 | + r += float64ToFloat32( ls * float64ToFloat32(hs + ss)); |
| 190 | + s2 = float64ToFloat32( hs * hs ); |
| 191 | + ht = float64ToFloat32( 3.0 + s2 + r ); |
| 192 | + ahs = fromWordf( ht ); |
| 193 | + ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ); |
| 194 | + ht = toWordf( ahs ); |
| 195 | + lt = float64ToFloat32( r - float64ToFloat32((ht-3.0) - s2)); |
| 196 | + |
| 197 | + // u+v = ss*(1+...): |
| 198 | + u = float64ToFloat32( hs * ht ); |
| 199 | + v = float64ToFloat32( float64ToFloat32( ls*ht ) + float64ToFloat32( lt*ss )); // eslint-disable-line max-len |
| 200 | + |
| 201 | + // 2/(3LN2) * (ss+...): |
| 202 | + hp = u + v; |
| 203 | + ahs = fromWordf( hp ); |
| 204 | + ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ); |
| 205 | + hp = toWordf( ahs ); |
| 206 | + lp = float64ToFloat32( v - float64ToFloat32(hp - u)); |
| 207 | + hz = float64ToFloat32( CP_HI * hp ); // CP_HI+CP_LO = 2/(3*LN2) |
| 208 | + lz = float64ToFloat32( float64ToFloat32( CP_LO*hp ) + float64ToFloat32( lp*CP ) + DP_LO[ k ]); // eslint-disable-line max-len |
| 209 | + |
| 210 | + // log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz |
| 211 | + dp = float64ToFloat32( DP_HI[ k ] ); |
| 212 | + t = float64ToFloat32( n ); |
| 213 | + // log2(ax) |
| 214 | + t1 = float64ToFloat32(( float64ToFloat32(hz+lz) + float64ToFloat32( dp )) + float64ToFloat32( t )); // eslint-disable-line max-len |
| 215 | + ahs = fromWordf( t1 ); |
| 216 | + ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ); |
| 217 | + t1 = toWordf( ahs ); |
| 218 | + t2 = float64ToFloat32( lz - float64ToFloat32(((t1-t) - dp) - float64ToFloat32( hz ))); // eslint-disable-line max-len |
| 219 | + |
| 220 | + return out; |
| 221 | +} |
| 222 | + |
| 223 | + |
| 224 | +// EXPORTS // |
| 225 | + |
| 226 | +module.exports = log2ax; |
0 commit comments