time-to-botec

log1p.js (9445B)

---

 /**
 * @license Apache-2.0
 *
 * Copyright (c) 2018 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.
 *
 *
 * ## Notice
 *
 * The following copyright, license, and long comment were part of the original implementation available as part of [FDLIBM]{@link http://www.netlib.org/fdlibm/s_log1p.c}. The implementation follows the original, but has been modified for JavaScript.
 *
 * ```text
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ```
 */
 
 'use strict';
 
 // MODULES //
 
 var isnan = require( './../../../../base/assert/is-nan' );
 var getHighWord = require( '@stdlib/number/float64/base/get-high-word' );
 var setHighWord = require( '@stdlib/number/float64/base/set-high-word' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 var FLOAT64_EXPONENT_BIAS = require( '@stdlib/constants/float64/exponent-bias' );
 var polyval = require( './polyval_lp.js' );
 
 
 // VARIABLES //
 
 // High and low words of ln(2):
 var LN2_HI = 6.93147180369123816490e-01; // 0x3fe62e42 0xfee00000
 var LN2_LO = 1.90821492927058770002e-10; // 0x3dea39ef 0x35793c76
 
 // sqrt(2)-1:
 var SQRT2M1 = 4.142135623730950488017e-01; // 0x3fda8279 0x99fcef34
 
 // sqrt(2)/2-1:
 var SQRT2HALFM1 = -2.928932188134524755992e-01; // 0xbfd2bec3 0x33018866
 
 // 2**-29:
 var SMALL = 1.862645149230957e-09; // 0x3e200000 0x00000000
 
 // 2**-54:
 var TINY = 5.551115123125783e-17;
 
 // Max integer (unsafe) => 2**53:
 var TWO53 = 9007199254740992;
 
 // 2/3:
 var TWO_THIRDS = 6.666666666666666666e-01;
 
 
 // MAIN //
 
 /**
 * Evaluates the natural logarithm of \\(1+x\\).
 *
 * ## Method
 *
 * 1.  Argument Reduction: find \\(k\\) and \\(f\\) such that
 *
 *     ```tex
 *     1+x = 2^k (1+f)
 *     ```
 *
 *     where
 *
 *     ```tex
 *     \frac{\sqrt{2}}{2} < 1+f < \sqrt{2}
 *     ```
 *
 *     <!-- <note> -->
 *
 *     If \\(k=0\\), then \\(f=x\\) is exact. However, if \\(k \neq 0\\), then \\(f\\) may not be representable exactly. In that case, a correction term is needed. Let
 *
 *     ```tex
 *     u = \operatorname{round}(1+x)
 *     ```
 *
 *     and
 *
 *     ```tex
 *     c = (1+x) - u
 *     ```
 *
 *     then
 *
 *     ```tex
 *     \ln (1+x) - \ln u \approx \frac{c}{u}
 *     ```
 *
 *     We can thus proceed to compute \\(\ln(u)\\), and add back the correction term \\(c/u\\).
 *
 *     <!-- </note> -->
 *
 *     <!-- <note> -->
 *
 *     When \\(x > 2^{53}\\), one can simply return \\(\ln(x)\\).
 *
 *     <!-- </note> -->
 *
 * 2.  Approximation of \\(\operatorname{log1p}(f)\\). Let
 *
 *     ```tex
 *     s = \frac{f}{2+f}
 *     ```
 *
 *     based on
 *
 *     ```tex
 *     \begin{align*}
 *     \ln 1+f &= \ln (1+s) - \ln (1-s) \\
 *             &= 2s + \frac{2}{3} s^3 + \frac{2}{5} s^5 + ... \\
 *             &= 2s + sR \\
 *     \end{align*}
 *     ```
 *
 *     We use a special Reme algorithm on \\(\[0,0.1716\]\\) to generate a polynomial of degree \\(14\\) to approximate \\(R\\). The maximum error of this polynomial approximation is bounded by \\(2^{-58.45}\\). In other words,
 *
 *     ```tex
 *     R(z) \approx \mathrm{Lp}_1 s^2 + \mathrm{Lp}_2 s^4 + \mathrm{Lp}_3 s^6 + \mathrm{Lp}_4 s^8 + \mathrm{Lp}_5 s^{10} + \mathrm{Lp}_6 s^{12} + \mathrm{Lp}_7 s^{14}
 *     ```
 *
 *     and
 *
 *     ```tex
 *     | \mathrm{Lp}_1 s^2 + \ldots + \mathrm{Lp}_7 s^14 - R(z) | \leq 2^{-58.45}
 *     ```
 *
 *     <!-- <note> -->
 *
 *     The values of \\(Lp1\\) to \\(Lp7\\) may be found in the source.
 *
 *     <!-- </note> -->
 *
 *     Note that
 *
 *     ```tex
 *     \begin{align*}
 *     2s &= f - sf \\
 *        &= f - \frac{f^2}{2} + s \frac{f^2}{2} \\
 *     \end{align*}
 *     ```
 *
 *     In order to guarantee error in \\(\ln\\) below \\(1\ \mathrm{ulp}\\), we compute the log by
 *
 *     ```tex
 *     \operatorname{log1p}(f) = f - \biggl(\frac{f^2}{2} - s\biggl(\frac{f^2}{2}+R\biggr)\biggr)
 *     ```
 *
 * 3.  Finally,
 *
 *     ```tex
 *     \begin{align*}
 *     \operatorname{log1p}(x) &= k \cdot \mathrm{ln2} + \operatorname{log1p}(f) \\
 *     &= k \cdot \mathrm{ln2}_{hi}+\biggl(f-\biggl(\frac{f^2}{2}-\biggl(s\biggl(\frac{f^2}{2}+R\biggr)+k \cdot \mathrm{ln2}_{lo}\biggr)\biggr)\biggr) \\
 *     \end{align*}
 *     ```
 *
 *     Here \\(\mathrm{ln2}\\) is split into two floating point numbers:
 *
 *     ```tex
 *     \mathrm{ln2}_{hi} + \mathrm{ln2}_{lo}
 *     ```
 *
 *     where \\(n \cdot \mathrm{ln2}_{hi}\\) is always exact for \\(|n| < 2000\\).
 *
 *
 * ## Special Cases
 *
 * -   \\(\operatorname{log1p}(x) = \mathrm{NaN}\\) with signal if \\(x < -1\\) (including \\(-\infty\\))
 * -   \\(\operatorname{log1p}(+\infty) = +\infty\\)
 * -   \\(\operatorname{log1p}(-1) = -\infty\\) with signal
 * -   \\(\operatorname{log1p}(\mathrm{NaN})= \mathrm{NaN}\\) with no signal
 *
 *
 * ## Notes
 *
 * -   According to an error analysis, the error is always less than \\(1\\) ulp (unit in the last place).
 *
 * -   The hexadecimal values are the intended ones for the used constants. The decimal values may be used, provided that the compiler will convert from decimal to binary accurately enough to produce the hexadecimal values shown.
 *
 * -   Assuming \\(\ln(x)\\) is accurate, the following algorithm can be used to evaluate \\(\operatorname{log1p}(x)\\) to within a few ULP:
 *
 *     ```javascript
 *     var u = 1.0 + x;
 *     if ( u === 1.0 ) {
 *         return x;
 *     } else {
 *         return ln(u) * (x/(u-1.0));
 *     }
 *     ```
 *
 *     See HP-15C Advanced Functions Handbook, p.193.
 *
 *
 * @param {number} x - input value
 * @returns {number} the natural logarithm of `1+x`
 *
 * @example
 * var v = log1p( 4.0 );
 * // returns ~1.609
 *
 * @example
 * var v = log1p( -1.0 );
 * // returns -Infinity
 *
 * @example
 * var v = log1p( 0.0 );
 * // returns 0.0
 *
 * @example
 * var v = log1p( -0.0 );
 * // returns -0.0
 *
 * @example
 * var v = log1p( -2.0 );
 * // returns NaN
 *
 * @example
 * var v = log1p( NaN );
 * // returns NaN
 */
 function log1p( x ) {
 	var hfsq;
 	var hu;
 	var y;
 	var f;
 	var c;
 	var s;
 	var z;
 	var R;
 	var u;
 	var k;
 
 	if ( x < -1.0 || isnan( x ) ) {
 		return NaN;
 	}
 	if ( x === -1.0 ) {
 		return NINF;
 	}
 	if ( x === PINF ) {
 		return x;
 	}
 	if ( x === 0.0 ) {
 		return x; // handle +-0 (IEEE 754-2008 spec)
 	}
 	// Set y = |x|:
 	if ( x < 0.0 ) {
 		y = -x;
 	} else {
 		y = x;
 	}
 	// Argument reduction...
 	k = 1;
 
 	// Check if argument reduction is needed and if we can just return a small value approximation requiring less computation but with equivalent accuracy...
 	if ( y < SQRT2M1 ) { // if |x| < sqrt(2)-1 => ~0.41422
 		if ( y < SMALL ) { // if |x| < 2**-29
 			if ( y < TINY ) { // if |x| < 2**-54
 				return x;
 			}
 			// Use a simple two-term Taylor series...
 			return x - ( x*x*0.5 );
 		}
 		// Check if `f=x` can be represented exactly (no need for correction terms), allowing us to bypass argument reduction...
 		if ( x > SQRT2HALFM1 ) { // if x > sqrt(2)/2-1 => ~-0.2929
 			// -0.2929 < x < 0.41422
 			k = 0;
 			f = x; // exact
 			hu = 1;
 		}
 	}
 	// Address case where `f` cannot be represented exactly...
 	if ( k !== 0 ) {
 		if ( y < TWO53 ) {
 			u = 1.0 + x;
 			hu = getHighWord( u );
 
 			// Bit shift to isolate the exponent and then subtract the bias:
 			k = (hu>>20) - FLOAT64_EXPONENT_BIAS;
 
 			// Correction term...
 			if ( k > 0 ) { // positive unbiased exponent
 				c = 1.0 - (u-x);
 			} else { // nonpositive unbiased exponent
 				c = x - (u-1.0);
 			}
 			c /= u;
 		} else {
 			u = x;
 			hu = getHighWord( u );
 
 			// Bit shift to isolate the exponent and then subtract the bias:
 			k = (hu>>20) - FLOAT64_EXPONENT_BIAS;
 
 			// Correction term is zero:
 			c = 0;
 		}
 		// Apply a bit mask (0 00000000000 11111111111111111111) to remove the exponent:
 		hu &= 0x000fffff; // max value => 1048575
 
 		// Check if u significand is less than sqrt(2) significand => 0x6a09e => 01101010000010011110
 		if ( hu < 434334 ) {
 			// Normalize u by setting the exponent to 1023 (bias) => 0x3ff00000 => 0 01111111111 00000000000000000000
 			u = setHighWord( u, hu|0x3ff00000 );
 		} else {
 			k += 1;
 
 			// Normalize u/2 by setting the exponent to 1022 (bias-1 => 2**-1 = 1/2) => 0x3fe00000 => 0 01111111110 00000000000000000000
 			u = setHighWord( u, hu|0x3fe00000 );
 
 			// Subtract hu significand from next largest hu => 0 00000000001 00000000000000000000 => 0x00100000 => 1048576
 			hu = (1048576-hu)>>2;
 		}
 		f = u - 1.0;
 	}
 	// Approximation of log1p(f)...
 	hfsq = 0.5 * f * f;
 	if ( hu === 0 ) { // if |f| < 2**-20
 		if ( f === 0.0 ) {
 			c += k * LN2_LO;
 			return ( k * LN2_HI ) + c;
 		}
 		R = hfsq * (1.0 - ( TWO_THIRDS*f ) ); // avoid division
 		return ( k*LN2_HI ) - ( (R - ( (k*LN2_LO) + c)) - f );
 	}
 	s = f / (2.0 + f);
 	z = s * s;
 
 	R = z * polyval( z );
 
 	if ( k === 0 ) {
 		return f - ( hfsq - ( s*(hfsq+R) ) );
 	}
 	return ( k*LN2_HI ) - ( (hfsq - ( (s*(hfsq+R)) + ((k*LN2_LO) + c))) - f );
 }
 
 
 // EXPORTS //
 
 module.exports = log1p;