|  | 
|  | 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 FLOAT32_EXPONENT_BIAS = require( '@stdlib/constants/float32/exponent-bias' ); | 
|  | 41 | +var FLOAT32_SIGNIFICAND_MASK = require( '@stdlib/constants/float32/significand-mask' ); | 
|  | 42 | +var polyvalL = require( '@stdlib/math/base/special/powf/lib/polyval_l.js' ); | 
|  | 43 | + | 
|  | 44 | + | 
|  | 45 | +// VARIABLES // | 
|  | 46 | + | 
|  | 47 | +// 0x0080000 = 524288 => 0 00000000000 10000000000000000000 | 
|  | 48 | +var HIGH_MIN_NORMAL_EXP = 0x00800000; // asm type annotation | 
|  | 49 | + | 
|  | 50 | +// 0x3f800000 = 1065353216 => 0 01111111000 00000000000000000000 | 
|  | 51 | +var HIGH_BIASED_EXP_0 = 0x3f800000; // asm type annotation | 
|  | 52 | + | 
|  | 53 | +// 0xfffff000 = 4294963200 => 1 11111111111 11111111111100000000 | 
|  | 54 | +var HIGH_BIASED_EXP_NEG_64 = 0xfffff000; // asm type annotation | 
|  | 55 | + | 
|  | 56 | +// 0x20000000 = 536870912 => 0 10000000000 00000000000000000000 | 
|  | 57 | +var HIGH_SIGNIFICAND_HALF = 0x20000000; // asm type annotation | 
|  | 58 | + | 
|  | 59 | +// 0x00400000 = 4194304 => 0 00000000000 10000000000000000000 | 
|  | 60 | +var HIGH_MIN_NORMAL_EXP_HALF = 0x00400000; // asm type annotation | 
|  | 61 | + | 
|  | 62 | +var HIGH_NUM_SIGNIFICAND_BITS = 23; // asm type annotation | 
|  | 63 | + | 
|  | 64 | +var TWO24 = 16777216.0;	// 0x4b800000 | 
|  | 65 | + | 
|  | 66 | +// 2/(3*LN2) | 
|  | 67 | +var CP = 9.6179670095e-01; // 0x3f76384f =2/(3ln2) | 
|  | 68 | + | 
|  | 69 | +// (float)CP | 
|  | 70 | +var CP_HI = 9.6191406250e-01; // 0x3f764000 =12b cp | 
|  | 71 | + | 
|  | 72 | +// Low: CP_HI | 
|  | 73 | +var CP_LO = -1.1736857402e-04; // 0xb8f623c6 =tail of CP_HI | 
|  | 74 | + | 
|  | 75 | +var BP = [ | 
|  | 76 | +	1.0, | 
|  | 77 | +	1.5 | 
|  | 78 | +]; | 
|  | 79 | +var DP_HI = [ | 
|  | 80 | +	0.0, | 
|  | 81 | +	5.84960938e-01 // 0x3f15c000 | 
|  | 82 | +]; | 
|  | 83 | +var DP_LO = [ | 
|  | 84 | +	0.0, | 
|  | 85 | +	1.56322085e-06 // 0x35d1cfdc | 
|  | 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 ahs; | 
|  | 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; // asm type annotation | 
|  | 129 | + | 
|  | 130 | +	// Check if `x` is subnormal... | 
|  | 131 | +	if ( ahx < HIGH_MIN_NORMAL_EXP ) { | 
|  | 132 | +		ax *= TWO24; | 
|  | 133 | +		n -= 24; // asm type annotation | 
|  | 134 | +		ahx = fromWordf( ax ); | 
|  | 135 | +	} | 
|  | 136 | +	// Extract the unbiased exponent of `x`: | 
|  | 137 | +	n += ((ahx >> HIGH_NUM_SIGNIFICAND_BITS) - FLOAT32_EXPONENT_BIAS); // asm type annotation | 
|  | 138 | + | 
|  | 139 | +	// Isolate the significand bits of `x`: | 
|  | 140 | +	j = (ahx & FLOAT32_SIGNIFICAND_MASK); // asm type annotation | 
|  | 141 | + | 
|  | 142 | +	// Normalize `ahx` by setting the (biased) exponent to `127`: | 
|  | 143 | +	ahx = (j | HIGH_BIASED_EXP_0); // asm type annotation | 
|  | 144 | + | 
|  | 145 | +	// Determine the interval of `|x|` by comparing significand bits... | 
|  | 146 | + | 
|  | 147 | +	// |x| < sqrt(3/2) | 
|  | 148 | +	if ( j <= 0x1cc471 ) { // 0 00000000001 11001100100010001110 | 
|  | 149 | +		k = 0; | 
|  | 150 | +	} | 
|  | 151 | +	// |x| < sqrt(3) | 
|  | 152 | +	else if ( j < 0x5db3d7 ) { // 0 00000000010 11101101100111110111 | 
|  | 153 | +		k = 1; | 
|  | 154 | +	} | 
|  | 155 | +	// |x| >= sqrtf(3) | 
|  | 156 | +	else { | 
|  | 157 | +		k = 0; | 
|  | 158 | +		n += 1; // asm type annotation | 
|  | 159 | +		ahx -= HIGH_MIN_NORMAL_EXP; | 
|  | 160 | +	} | 
|  | 161 | +	// Load the normalized float word into `|x|`: | 
|  | 162 | +	ahx = float64ToFloat32( ahx ); | 
|  | 163 | +	ax = toWordf( ahx ); | 
|  | 164 | + | 
|  | 165 | +	// Compute `ss = hs + ls = (x-1)/(x+1)` or `(x-1.5)/(x+1.5)`: | 
|  | 166 | +	bp = float64ToFloat32( BP[ k ] ); // BP[0] = 1.0, BP[1] = 1.5 | 
|  | 167 | +	u = float64ToFloat32( ax - bp ); // (x-1) || (x-1.5) | 
|  | 168 | +	v = float64ToFloat32( 1.0 / float64ToFloat32(ax + bp) ); // 1/(x+1) || 1/(x+1.5) | 
|  | 169 | +	ss = float64ToFloat32( u * v ); | 
|  | 170 | +	hs = float64ToFloat32( ss ); | 
|  | 171 | +	ahs = fromWordf( hs ); | 
|  | 172 | +	ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ) | 
|  | 173 | +	hs = toWordf( ahs ); | 
|  | 174 | + | 
|  | 175 | +	// Compute `ht = ax + bp` High | 
|  | 176 | +	ahs = ((ahx>>1) & HIGH_BIASED_EXP_NEG_64) | HIGH_SIGNIFICAND_HALF; | 
|  | 177 | +	ahs += HIGH_MIN_NORMAL_EXP_HALF+(k << 21); // `(k<<21)` can be considered the word equivalent of `1.0` or `1.5` | 
|  | 178 | +	ahs = float64ToFloat32( ahs ); | 
|  | 179 | +	ht = toWordf( ahs ); | 
|  | 180 | +	lt = float64ToFloat32( ax - float64ToFloat32(ht - bp) ); | 
|  | 181 | +	ls = float64ToFloat32( v * ( float64ToFloat32( u - float64ToFloat32(hs*ht) ) - float64ToFloat32( hs*lt ) ) ); // eslint-disable-line max-len | 
|  | 182 | + | 
|  | 183 | +	// Compute `log(ax)`... | 
|  | 184 | + | 
|  | 185 | +	s2 = float64ToFloat32( ss * ss); | 
|  | 186 | +	r = float64ToFloat32( s2 * s2 * float64ToFloat32( polyvalL( s2 ))); | 
|  | 187 | +	r += float64ToFloat32( ls * float64ToFloat32(hs + ss)); | 
|  | 188 | +	s2 = float64ToFloat32( hs * hs ); | 
|  | 189 | +	ht = float64ToFloat32( 3.0 + s2 + r ); | 
|  | 190 | +	ahs = fromWordf( ht ); | 
|  | 191 | +	ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ); | 
|  | 192 | +	ht = toWordf( ahs ); | 
|  | 193 | +	lt = float64ToFloat32( r - float64ToFloat32((ht-3.0) - s2)); | 
|  | 194 | + | 
|  | 195 | +	// u+v = ss*(1+...): | 
|  | 196 | +	u = float64ToFloat32( hs * ht ); | 
|  | 197 | +	v = float64ToFloat32( float64ToFloat32( ls*ht ) + float64ToFloat32( lt*ss )); // eslint-disable-line max-len | 
|  | 198 | + | 
|  | 199 | +	// 2/(3LN2) * (ss+...): | 
|  | 200 | +	hp = u + v; | 
|  | 201 | +	ahs = fromWordf( hp ); | 
|  | 202 | +	ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ); | 
|  | 203 | +	hp = toWordf( ahs ); | 
|  | 204 | +	lp = float64ToFloat32( v - float64ToFloat32(hp - u)); | 
|  | 205 | +	hz = float64ToFloat32( CP_HI * hp ); // CP_HI+CP_LO = 2/(3*LN2) | 
|  | 206 | +	lz = float64ToFloat32( float64ToFloat32( CP_LO*hp ) + float64ToFloat32( lp*CP ) + DP_LO[ k ]); // eslint-disable-line max-len | 
|  | 207 | + | 
|  | 208 | +	// log2(ax) = (ss+...)*2/(3*LN2) = n + dp + hz + lz | 
|  | 209 | +	dp = float64ToFloat32( DP_HI[ k ] ); | 
|  | 210 | +	t = float64ToFloat32( n ); | 
|  | 211 | +	// log2(ax) | 
|  | 212 | +	t1 = float64ToFloat32(( float64ToFloat32(hz+lz) + float64ToFloat32( dp )) + float64ToFloat32( t )); // eslint-disable-line max-len | 
|  | 213 | +	ahs = fromWordf( t1 ); | 
|  | 214 | +	ahs = float64ToFloat32( ahs & HIGH_BIASED_EXP_NEG_64 ); | 
|  | 215 | +	t1 = toWordf( ahs ); | 
|  | 216 | +	t2 = float64ToFloat32( lz - float64ToFloat32(((t1-t) - dp) - float64ToFloat32( hz ))); // eslint-disable-line max-len | 
|  | 217 | + | 
|  | 218 | +	return out; | 
|  | 219 | +} | 
|  | 220 | + | 
|  | 221 | + | 
|  | 222 | +// EXPORTS // | 
|  | 223 | + | 
|  | 224 | +module.exports = log2ax; | 
0 commit comments