feat: port exp2 from OpenLibm

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 <das@FreeBSD.ORG>
+* 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
 }