Skip to content

Commit be5ece8

Browse files
gunjjoshikgryte
andauthored
refactor: update math/base/special/log2 implementation
PR-URL: #2172 Closes: #2170 Co-authored-by: Athan Reines <[email protected]> Reviewed-by: Athan Reines <[email protected]>
1 parent 14b3434 commit be5ece8

File tree

11 files changed

+98
-498
lines changed

11 files changed

+98
-498
lines changed

lib/node_modules/@stdlib/math/base/special/log2/lib/klog.js

-102
This file was deleted.

lib/node_modules/@stdlib/math/base/special/log2/lib/main.js

+40-10
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
*
1919
* ## Notice
2020
*
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_log2.c}. The implementation follows the original, but has been modified for JavaScript.
21+
* The following copyright and license were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/12.2.0/lib/msun/src/e_log2.c}. The implementation follows the original, but has been modified for JavaScript.
2222
*
2323
* ```text
2424
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -43,7 +43,7 @@ var ABS_MASK = require( '@stdlib/constants/float64/high-word-abs-mask' );
4343
var HIGH_SIGNIFICAND_MASK = require( '@stdlib/constants/float64/high-word-significand-mask' );
4444
var BIAS = require( '@stdlib/constants/float64/exponent-bias' );
4545
var NINF = require( '@stdlib/constants/float64/ninf' );
46-
var klog = require( './klog.js' );
46+
var kernelLog1p = require( '@stdlib/math/base/special/kernel-log1p' );
4747

4848

4949
// VARIABLES //
@@ -98,11 +98,17 @@ var WORDS = [ 0|0, 0|0 ];
9898
* // returns NaN
9999
*/
100100
function log2( x ) {
101-
var hi;
102-
var lo;
101+
var valHi;
102+
var valLo;
103+
var hfsq;
103104
var hx;
104105
var lx;
106+
var hi;
107+
var lo;
105108
var f;
109+
var R;
110+
var w;
111+
var y;
106112
var i;
107113
var k;
108114

@@ -127,18 +133,42 @@ function log2( x ) {
127133
if ( hx >= HIGH_MAX_NORMAL_EXP ) {
128134
return x + x;
129135
}
136+
// Case: log(1) = +0
137+
if ( hx === HIGH_BIASED_EXP_0 && lx === 0 ) {
138+
return 0.0;
139+
}
130140
k += ( (hx>>20) - BIAS )|0; // asm type annotation
131141
hx &= HIGH_SIGNIFICAND_MASK;
132-
i = ( ( hx+0x95f64 ) & 0x100000 )|0; // asm type annotation
142+
i = ( ( hx+0x95f64 ) & HIGH_MIN_NORMAL_EXP )|0; // asm type annotation
133143

134144
// Normalize x or x/2...
135145
x = setHighWord( x, hx|(i^HIGH_BIASED_EXP_0) );
136146
k += (i>>20)|0; // asm type annotation
137-
f = klog( x );
138-
x -= 1;
139-
hi = setLowWord( x, 0 );
140-
lo = x - hi;
141-
return ( (x+f)*IVLN2LO ) + ( (lo+f)*IVLN2HI ) + ( hi*IVLN2HI ) + k;
147+
y = k;
148+
f = x - 1.0;
149+
hfsq = 0.5 * f * f;
150+
R = kernelLog1p( f );
151+
152+
/*
153+
* Notes:
154+
*
155+
* - `f-hfsq` must (for args near `1`) be evaluated in extra precision to avoid a large cancellation when `x` is near `sqrt(2)` or `1/sqrt(2)`.This is fairly efficient since `f-hfsq` only depends on `f`, so can be evaluated in parallel with `R`. Not combining `hfsq` with `R` also keeps `R` small (though not as small as a true `lo` term would be), so that extra precision is not needed for terms involving `R`.
156+
* - When implemented in C, compiler bugs involving extra precision used to break Dekker's theorem for spitting `f-hfsq` as `hi+lo`. These problems are now automatically avoided as a side effect of the optimization of combining the Dekker splitting step with the clear-low-bits step.
157+
* - `y` must (for args near `sqrt(2)` and `1/sqrt(2)`) be added in extra precision to avoid a very large cancellation when `x` is very near these values. Unlike the above cancellations, this problem is specific to base `2`. It is strange that adding `+-1` is so much harder than adding `+-ln2` or `+-log10_2`.
158+
* - This implementation uses Dekker's theorem to normalize `y+val_hi`, so, when implemented in C, compiler bugs may reappear in some configurations.
159+
* - The multi-precision calculations for the multiplications are routine.
160+
*/
161+
hi = f - hfsq;
162+
hi = setLowWord( hi, 0 );
163+
lo = ( f - hi ) - hfsq + R;
164+
valHi = hi * IVLN2HI;
165+
valLo = ( ( lo + hi ) * IVLN2LO ) + ( lo * IVLN2HI );
166+
167+
w = y + valHi;
168+
valLo += ( y - w ) + valHi;
169+
valHi = w;
170+
171+
return valLo + valHi;
142172
}
143173

144174

lib/node_modules/@stdlib/math/base/special/log2/lib/polyval_p.js

-47
This file was deleted.

lib/node_modules/@stdlib/math/base/special/log2/lib/polyval_q.js

-47
This file was deleted.

lib/node_modules/@stdlib/math/base/special/log2/manifest.json

+6-3
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,8 @@
4545
"@stdlib/constants/float64/high-word-abs-mask",
4646
"@stdlib/constants/float64/high-word-significand-mask",
4747
"@stdlib/constants/float64/exponent-bias",
48-
"@stdlib/constants/float64/ninf"
48+
"@stdlib/constants/float64/ninf",
49+
"@stdlib/math/base/special/kernel-log1p"
4950
]
5051
},
5152
{
@@ -67,7 +68,8 @@
6768
"@stdlib/constants/float64/high-word-abs-mask",
6869
"@stdlib/constants/float64/high-word-significand-mask",
6970
"@stdlib/constants/float64/exponent-bias",
70-
"@stdlib/constants/float64/ninf"
71+
"@stdlib/constants/float64/ninf",
72+
"@stdlib/math/base/special/kernel-log1p"
7173
]
7274
},
7375
{
@@ -89,7 +91,8 @@
8991
"@stdlib/constants/float64/high-word-abs-mask",
9092
"@stdlib/constants/float64/high-word-significand-mask",
9193
"@stdlib/constants/float64/exponent-bias",
92-
"@stdlib/constants/float64/ninf"
94+
"@stdlib/constants/float64/ninf",
95+
"@stdlib/math/base/special/kernel-log1p"
9396
]
9497
}
9598
]

lib/node_modules/@stdlib/math/base/special/log2/package.json

-1
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,6 @@
2121
"example": "./examples",
2222
"include": "./include",
2323
"lib": "./lib",
24-
"scripts": "./scripts",
2524
"src": "./src",
2625
"test": "./test"
2726
},

0 commit comments

Comments
 (0)