From 37cc5467b271641c8dc5e4e8b42737f78b15e327 Mon Sep 17 00:00:00 2001 From: Kamal Singh Rautela Date: Tue, 10 Feb 2026 11:59:04 +0530 Subject: [PATCH] fix: initialize hx in exp2 and refactor implementation --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../math/base/special/exp2/lib/main.js | 322 ++++++++++++++---- .../math/base/special/exp2/lib/polyval_p.js | 47 --- .../math/base/special/exp2/lib/polyval_q.js | 47 --- 3 files changed, 262 insertions(+), 154 deletions(-) delete mode 100644 lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_p.js delete mode 100644 lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_q.js diff --git a/lib/node_modules/@stdlib/math/base/special/exp2/lib/main.js b/lib/node_modules/@stdlib/math/base/special/exp2/lib/main.js index ed01a638747f..e25ce8ce9e7c 100644 --- a/lib/node_modules/@stdlib/math/base/special/exp2/lib/main.js +++ b/lib/node_modules/@stdlib/math/base/special/exp2/lib/main.js @@ -18,62 +18,228 @@ * * ## Notice * -* The original C code, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. The implementation follows the original, but has been modified for JavaScript. +* The following copyright, license, and long comment were part of the original implementation available as part of [FreeBSD]{@link https://svnweb.freebsd.org/base/release/12.2.0/lib/msun/src/s_exp2.c}. The implementation follows the original, but has been modified for JavaScript. * * ```text -* Copyright 1984, 1995, 2000 by Stephen L. Moshier +* Copyright (c) 2005 David Schultz +* All rights reserved. * -* Some software in this archive may be from the book _Methods and Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster International, 1989) or from the Cephes Mathematical Library, a commercial product. In either event, it is copyrighted by the author. What you see here may be used freely but it comes with no support or guarantee. +* Redistribution and use in source and binary forms, with or without +* modification, are permitted provided that the following conditions +* are met: +* 1. Redistributions of source code must retain the above copyright +* notice, this list of conditions and the following disclaimer. +* 2. Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions and the following disclaimer in the +* documentation and/or other materials provided with the distribution. * -* Stephen L. Moshier -* moshier@na-net.ornl.gov +* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND +* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE +* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +* SUCH DAMAGE. * ``` */ -'use strict'; +/* + * exp2(x): compute the base 2 exponential of x + * + * Accuracy: Peak error < 0.503 ulp for normalized results. + * + * Method: (accurate tables) + * + * Reduce x: + * x = k + y, for integer k and |y| <= 1/2. + * Thus we have exp2(x) = 2**k * exp2(y). + * + * Reduce y: + * y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE. + * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]), + * with |z - eps[i]| <= 2**-9 + 2**-39 for the table used. + * + * We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via + * a degree-5 minimax polynomial with maximum error under 1.3 * 2**-61. + * The values in exp2t[] and eps[] are chosen such that + * exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such + * that exp2t[i] is accurate to 2**-64. + * + * Note that the range of i is +-TBLSIZE/2, so we actually index the tables + * by i0 = i + TBLSIZE/2. For cache efficiency, exp2t[] and eps[] are + * virtual tables, interleaved in the real table tbl[]. + * + * This method is due to Gal, with many details due to Gal and Bachelis: + * + * Gal, S. and Bachelis, B. An Accurate Elementary Mathematical Library + * for the IEEE Floating Point Standard. TOMS 17(1), 26-46 (1991). + */ -// TODO: replace with TOMS (Openlibm) algo (updating license header and long comment) +'use strict'; // MODULES // -var FLOAT64_MAX_BASE2_EXPONENT = require( '@stdlib/constants/float64/max-base2-exponent' ); // eslint-disable-line id-length -var FLOAT64_MIN_BASE2_EXPONENT = require( '@stdlib/constants/float64/min-base2-exponent' ); // eslint-disable-line id-length -var round = require( '@stdlib/math/base/special/round' ); -var ldexp = require( '@stdlib/math/base/special/ldexp' ); -var isnan = require( '@stdlib/math/base/assert/is-nan' ); -var PINF = require( '@stdlib/constants/float64/pinf' ); -var polyvalP = require( './polyval_p.js' ); -var polyvalQ = require( './polyval_q.js' ); +var setHighWord = require( '@stdlib/number/float64/base/set-high-word' ); +var getHighWord = require( '@stdlib/number/float64/base/get-high-word' ); +var getLowWord = require( '@stdlib/number/float64/base/get-low-word' ); +var BIAS = require( '@stdlib/constants/float64/exponent-bias' ); + + +// VARIABLES // + +// TBLSIZE = 1 << TBLBITS where TBLBITS = 6 +var TBLSIZE = 64; +var TBLBITS = 6; +var REDUX = 6755399441055744.0; // 0x1.8p52 +var P1 = 0.6931471805599453; // 0x1.62e42fefa39efp-1 +var P2 = 0.2402265069591007; // 0x1.ebfbdff82c58fp-3 +var P3 = 0.05550410866482158; // 0x1.c6b08d704a0bfp-5 +var P4 = 0.009618129107628477; // 0x1.3b2ab6fba4e77p-7 +var P5 = 0.0013333558146428445; // 0x1.5d87fe781d948p-10 + +// Table of values for 2^(i/TBLSIZE) interleaved with small correction terms eps[i]. + +// Format: { expected result, epsilon correction } +var TBL = [ + 1.00000000000000000e+00, + 0.00000000000000000e+00, + 1.01088928605170048e+00, + 2.08166817117216851e-17, + 1.02189714865411663e+00, + -7.28583859910258980e-17, + 1.03302487902122841e+00, + -1.38777878078144568e-17, + 1.04427378242741375e+00, + -1.17961196366422882e-16, + 1.05564517836055716e+00, + 0.00000000000000000e+00, + 1.06714040067682370e+00, + 1.11022302462515654e-16, + 1.07876079775711986e+00, + 8.32667268468867405e-17, + 1.09050773266525769e+00, + 2.77555756156289135e-17, + 1.10238258330784089e+00, + -5.55111512312578270e-17, + 1.11438674259589243e+00, + -1.38777878078144568e-16, + 1.12652161860824185e+00, + -5.55111512312578270e-17, + 1.13878863475669156e+00, + -1.11022302462515654e-16, + 1.15118922995298267e+00, + -2.77555756156289135e-17, + 1.16372485877757748e+00, + -5.55111512312578270e-17, + 1.17639699165028122e+00, + -5.55111512312578270e-17, + 1.18920711500272103e+00, + -5.55111512312578270e-17, + 1.20215673145270308e+00, + -5.55111512312578270e-17, + 1.21524735998046896e+00, + 1.11022302462515654e-16, + 1.22848053610687002e+00, + 0.00000000000000000e+00, + 1.24185781207348400e+00, + -5.55111512312578270e-17, + 1.25538075702469110e+00, + 0.00000000000000000e+00, + 1.26905095719173322e+00, + 0.00000000000000000e+00, + 1.28287001607877826e+00, + 0.00000000000000000e+00, + 1.29683955465100964e+00, + -5.55111512312578270e-17, + 1.31096121152476441e+00, + 5.55111512312578270e-17, + 1.32523664315974132e+00, + 5.55111512312578270e-17, + 1.33966752405330292e+00, + -1.11022302462515654e-16, + 1.35425554693689265e+00, + -5.55111512312578270e-17, + 1.36900242297459052e+00, + -1.11022302462515654e-16, + 1.38390988196383202e+00, + 5.55111512312578270e-17, + 1.39897967253831124e+00, + 1.11022302462515654e-16, + 1.41421356237309515e+00, + 1.11022302462515654e-16, + 1.42961333839197002e+00, + 0.00000000000000000e+00, + 1.44518080697704665e+00, + 0.00000000000000000e+00, + 1.46091779418064704e+00, + 0.00000000000000000e+00, + 1.47682614593949935e+00, + 0.00000000000000000e+00, + 1.49290772829126484e+00, + 0.00000000000000000e+00, + 1.50916442759342284e+00, + 1.11022302462515654e-16, + 1.52559815074453842e+00, + 1.11022302462515654e-16, + 1.54221082540794074e+00, + -1.11022302462515654e-16, + 1.55900440023783693e+00, + 0.00000000000000000e+00, + 1.57598084510788650e+00, + 0.00000000000000000e+00, + 1.59314215134226700e+00, + 1.11022302462515654e-16, + 1.61049033194925428e+00, + 0.00000000000000000e+00, + 1.62802742185734783e+00, + 1.11022302462515654e-16, + 1.64575547815396495e+00, + 1.11022302462515654e-16, + 1.66367658032673638e+00, + 0.00000000000000000e+00, + 1.68179283050742900e+00, + -1.11022302462515654e-16, + 1.70010635371852348e+00, + 0.00000000000000000e+00, + 1.71861929812247793e+00, + 0.00000000000000000e+00, + 1.73733383527370622e+00, + 0.00000000000000000e+00, + 1.75625216037329945e+00, + 0.00000000000000000e+00, + 1.77537649252652119e+00, + 0.00000000000000000e+00, + 1.79470907500310717e+00, + 0.00000000000000000e+00, + 1.81425217550039886e+00, + 1.11022302462515654e-16, + 1.83400808640934243e+00, + 0.00000000000000000e+00, + 1.85397912508338547e+00, + -1.11022302462515654e-16, + 1.87416763411029996e+00, + 0.00000000000000000e+00, + 1.89457598158696561e+00, + 0.00000000000000000e+00, + 1.91520656139714740e+00, + 1.11022302462515654e-16, + 1.93606179349229435e+00, + -1.11022302462515654e-16, + 1.95714412417540018e+00, + -1.11022302462515654e-16, + 1.97845602638795093e+00, + 0.00000000000000000e+00 +]; // MAIN // /** -* Evaluates the base `2` exponential function. -* -* ## Method -* -* - Range reduction is accomplished by separating the argument into an integer \\( k \\) and fraction \\( f \\) such that -* -* ```tex -* 2^x = 2^k 2^f -* ``` -* -* - A Pade' approximate -* -* ```tex -* 1 + 2x \frac{\mathrm{P}\left(x^2\right)}{\mathrm{Q}\left(x^2\right) - x \mathrm{P}\left(x^2\right)} -* ``` -* -* approximates \\( 2^x \\) in the basic range \\( \[-0.5, 0.5] \\). -* -* ## Notes -* -* - Relative error: -* -* | arithmetic | domain | # trials | peak | rms | -* |:----------:|:-----------:|:--------:|:-------:|:-------:| -* | IEEE | -1022,+1024 | 30000 | 1.8e-16 | 5.4e-17 | +* Evaluates the base 2 exponential function. * * @param {number} x - input value * @returns {number} function value @@ -94,30 +260,66 @@ var polyvalQ = require( './polyval_q.js' ); * var v = exp2( NaN ); * // returns NaN */ -function exp2( x ) { - var px; - var xx; - var n; - if ( isnan( x ) ) { - return x; - } - if ( x > FLOAT64_MAX_BASE2_EXPONENT ) { - return PINF; +function exp2(x) { + var twopk; + var hx = getHighWord(x); + var lx; + var i0; + var ix = hx & 0x7fffffff; + var k; + var t; + var z; + var r; + + if (ix >= 0x40900000) { // |x| >= 1024 + if (ix >= 0x7ff00000) { // |x| >= Inf + lx = getLowWord(x); + if (((ix & 0xfffff) | lx) !== 0 || (hx & 0x80000000) === 0) { + return x + x; // NaN or +Inf + } + return 0.0; // -Inf + } + if (x >= 1024.0) { + return Number.POSITIVE_INFINITY; // overflow + } + if (x <= -1024.0) { + return 0.0; // underflow + } + } else if (ix < 0x3c900000) { // |x| < 2^-54 + return 1.0 + x; } - if ( x < FLOAT64_MIN_BASE2_EXPONENT ) { - return 0.0; + // Reduce x, computing z, i0, and k... + t = (x * TBLSIZE) + REDUX; + i0 = getLowWord(t); + k = (i0 >> TBLBITS) << 20; + i0 = (i0 & (TBLSIZE - 1)) << 1; + t -= REDUX; + t /= TBLSIZE; + z = x - t; + + // Compute r = exp2(y) = exp2t[i0] * p(z - eps[i])... + t = TBL[i0]; // exp2t[i0] + z -= TBL[i0 + 1]; // eps[i0] + + if (k >= -1070596096) { // k >= -(1021 << 20) + twopk = setHighWord(0.0, (BIAS << 20) + k); + } else { + twopk = setHighWord(0.0, ((BIAS + 1000) << 20) + k); } - // Separate into integer and fractional parts... - n = round( x ); - x -= n; - xx = x * x; - px = x * polyvalP( xx ); - x = px / ( polyvalQ( xx ) - px ); - x = 1.0 + ldexp( x, 1 ); + // Pade approximation... + r = t + (t * z * (P1 + (z * (P2 + (z * (P3 + (z * (P4 + (z * P5))))))))); - // Scale by power of 2: - return ldexp( x, n ); + // Scale by 2**(k>>20)... + if (k >= -1070596096) { // k >= -(1021 << 20) + if (k === 1073741824) { // k == 1024 << 20 + return r * 2.0 * 8.98846567431158e307; // 0x1p1023 + } + // eslint-disable-next-line @cspell/spellchecker + return r * twopk; // twopk is 2^k + } + // eslint-disable-next-line @cspell/spellchecker + return r * twopk * 9.332636185032189e-302; // twom1000 } diff --git a/lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_p.js b/lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_p.js deleted file mode 100644 index 5c7eb002067d..000000000000 --- a/lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_p.js +++ /dev/null @@ -1,47 +0,0 @@ -/** -* @license Apache-2.0 -* -* Copyright (c) 2024 The Stdlib Authors. -* -* Licensed under the Apache License, Version 2.0 (the "License"); -* you may not use this file except in compliance with the License. -* You may obtain a copy of the License at -* -* http://www.apache.org/licenses/LICENSE-2.0 -* -* Unless required by applicable law or agreed to in writing, software -* distributed under the License is distributed on an "AS IS" BASIS, -* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -* See the License for the specific language governing permissions and -* limitations under the License. -*/ - -/* This is a generated file. Do not edit directly. */ -'use strict'; - -// MAIN // - -/** -* Evaluates a polynomial. -* -* ## Notes -* -* - The implementation uses [Horner's rule][horners-method] for efficient computation. -* -* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method -* -* @private -* @param {number} x - value at which to evaluate the polynomial -* @returns {number} evaluated polynomial -*/ -function evalpoly( x ) { - if ( x === 0.0 ) { - return 1513.906801156151; - } - return 1513.906801156151 + (x * (20.202065669316532 + (x * 0.023093347705734523))); // eslint-disable-line max-len -} - - -// EXPORTS // - -module.exports = evalpoly; diff --git a/lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_q.js b/lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_q.js deleted file mode 100644 index 69f4d97f7f1f..000000000000 --- a/lib/node_modules/@stdlib/math/base/special/exp2/lib/polyval_q.js +++ /dev/null @@ -1,47 +0,0 @@ -/** -* @license Apache-2.0 -* -* Copyright (c) 2024 The Stdlib Authors. -* -* Licensed under the Apache License, Version 2.0 (the "License"); -* you may not use this file except in compliance with the License. -* You may obtain a copy of the License at -* -* http://www.apache.org/licenses/LICENSE-2.0 -* -* Unless required by applicable law or agreed to in writing, software -* distributed under the License is distributed on an "AS IS" BASIS, -* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -* See the License for the specific language governing permissions and -* limitations under the License. -*/ - -/* This is a generated file. Do not edit directly. */ -'use strict'; - -// MAIN // - -/** -* Evaluates a polynomial. -* -* ## Notes -* -* - The implementation uses [Horner's rule][horners-method] for efficient computation. -* -* [horners-method]: https://en.wikipedia.org/wiki/Horner%27s_method -* -* @private -* @param {number} x - value at which to evaluate the polynomial -* @returns {number} evaluated polynomial -*/ -function evalpoly( x ) { - if ( x === 0.0 ) { - return 4368.211668792106; - } - return 4368.211668792106 + (x * (233.1842117223149 + (x * 1.0))); -} - - -// EXPORTS // - -module.exports = evalpoly;