time-to-botec

pow.js (10177B)

---

 /**
 * @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 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/s_pow.c}. The implementation follows the original, but has been modified for JavaScript.
 *
 * ```text
 * Copyright (C) 2004 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 isOdd = require( './../../../../base/assert/is-odd' );
 var isInfinite = require( './../../../../base/assert/is-infinite' );
 var isInteger = require( './../../../../base/assert/is-integer' );
 var sqrt = require( './../../../../base/special/sqrt' );
 var abs = require( './../../../../base/special/abs' );
 var toWords = require( '@stdlib/number/float64/base/to-words' );
 var setLowWord = require( '@stdlib/number/float64/base/set-low-word' );
 var uint32ToInt32 = require( '@stdlib/number/uint32/base/to-int32' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var xIsZero = require( './x_is_zero.js' );
 var yIsHuge = require( './y_is_huge.js' );
 var yIsInfinite = require( './y_is_infinite.js' );
 var log2ax = require( './log2ax.js' );
 var logx = require( './logx.js' );
 var pow2 = require( './pow2.js' );
 
 
 // VARIABLES //
 
 // 0x7fffffff = 2147483647 => 0 11111111111 11111111111111111111
 var ABS_MASK = 0x7fffffff|0; // asm type annotation
 
 // 0x3fefffff = 1072693247 => 0 01111111110 11111111111111111111 => biased exponent: 1022 = -1+1023 => 2^-1
 var HIGH_MAX_NEAR_UNITY = 0x3fefffff|0; // asm type annotation
 
 // 0x41e00000 = 1105199104 => 0 10000011110 00000000000000000000 => biased exponent: 1054 = 31+1023 => 2^31
 var HIGH_BIASED_EXP_31 = 0x41e00000|0; // asm type annotation
 
 // 0x43f00000 = 1139802112 => 0 10000111111 00000000000000000000 => biased exponent: 1087 = 64+1023 => 2^64
 var HIGH_BIASED_EXP_64 = 0x43f00000|0; // asm type annotation
 
 // 0x40900000 = 1083179008 => 0 10000001001 00000000000000000000 => biased exponent: 1033 = 10+1023 => 2^10 = 1024
 var HIGH_BIASED_EXP_10 = 0x40900000|0; // asm type annotation
 
 // 0x3ff00000 = 1072693248 => 0 01111111111 00000000000000000000 => biased exponent: 1023 = 0+1023 => 2^0 = 1
 var HIGH_BIASED_EXP_0 = 0x3ff00000|0; // asm type annotation
 
 // 0x4090cc00 = 1083231232 => 0 10000001001 00001100110000000000
 var HIGH_1075 = 0x4090cc00|0; // asm type annotation
 
 // 0xc090cc00 = 3230714880 => 1 10000001001 00001100110000000000
 var HIGH_NEG_1075 = 0xc090cc00>>>0; // asm type annotation
 
 var HIGH_NUM_NONSIGN_BITS = 31|0; // asm type annotation
 
 var HUGE = 1.0e300;
 var TINY = 1.0e-300;
 
 // -(1024-log2(ovfl+.5ulp))
 var OVT = 8.0085662595372944372e-17;
 
 // High/low words workspace:
 var WORDS = [ 0|0, 0|0 ]; // WARNING: not thread safe
 
 // Log workspace:
 var LOG_WORKSPACE = [ 0.0, 0.0 ]; // WARNING: not thread safe
 
 
 // MAIN //
 
 /**
 * Evaluates the exponential function.
 *
 * ## Method
 *
 * 1.  Let \\(x = 2^n (1+f)\\).
 *
 * 2.  Compute \\(\operatorname{log2}(x)\\) as
 *
 *     ```tex
 *     \operatorname{log2}(x) = w_1 + w_2
 *     ```
 *
 *     where \\(w_1\\) has \\(53 - 24 = 29\\) bit trailing zeros.
 *
 * 3.  Compute
 *
 *     ```tex
 *     y \cdot \operatorname{log2}(x) = n + y^\prime
 *     ```
 *
 *     by simulating multi-precision arithmetic, where \\(|y^\prime| \leq 0.5\\).
 *
 * 4.  Return
 *
 *     ```tex
 *     x^y = 2^n e^{y^\prime \cdot \mathrm{log2}}
 *     ```
 *
 * ## Special Cases
 *
 * ```tex
 * \begin{align*}
 * x^{\mathrm{NaN}} &= \mathrm{NaN} & \\
 * (\mathrm{NaN})^y &= \mathrm{NaN} & \\
 * 1^y &= 1 & \\
 * x^0 &= 1 & \\
 * x^1 &= x & \\
 * (\pm 0)^\infty &= +0 & \\
 * (\pm 0)^{-\infty} &= +\infty & \\
 * (+0)^y &= +0 & \mathrm{if}\ y > 0 \\
 * (+0)^y &= +\infty & \mathrm{if}\ y < 0 \\
 * (-0)^y &= -\infty & \mathrm{if}\ y\ \mathrm{is\ an\ odd\ integer\ and}\ y < 0 \\
 * (-0)^y &= +\infty & \mathrm{if}\ y\ \mathrm{is\ not\ an\ odd\ integer\ and}\ y < 0 \\
 * (-0)^y &= -0 & \mathrm{if}\ y\ \mathrm{is\ an\ odd\ integer\ and}\ y > 0 \\
 * (-0)^y &= +0 & \mathrm{if}\ y\ \mathrm{is\ not\ an\ odd\ integer\ and}\ y > 0 \\
 * (-1)^{\pm\infty} &= \mathrm{NaN} & \\
 * x^{\infty} &= +\infty & |x| > 1 \\
 * x^{\infty} &= +0 & |x| < 1 \\
 * x^{-\infty} &= +0 & |x| > 1 \\
 * x^{-\infty} &= +\infty & |x| < 1 \\
 * (-\infty)^y &= (-0)^y & \\
 * \infty^y &= +0 & y < 0 \\
 * \infty^y &= +\infty & y > 0 \\
 * x^y &= \mathrm{NaN} & \mathrm{if}\ y\ \mathrm{is\ not\ a\ finite\ integer\ and}\ x < 0
 * \end{align*}
 * ```
 *
 * ## Notes
 *
 * -   \\(\operatorname{pow}(x,y)\\) returns \\(x^y\\) nearly rounded. In particular, \\(\operatorname{pow}(<\mathrm{integer}>,<\mathrm{integer}>)\\) **always** returns the correct integer, provided the value is representable.
 * -   The hexadecimal values shown in the source code are the intended values for used constants. Decimal values may be used, provided the compiler will accurately convert decimal to binary in order to produce the hexadecimal values.
 *
 *
 * @param {number} x - base
 * @param {number} y - exponent
 * @returns {number} function value
 *
 * @example
 * var v = pow( 2.0, 3.0 );
 * // returns 8.0
 *
 * @example
 * var v = pow( 4.0, 0.5 );
 * // returns 2.0
 *
 * @example
 * var v = pow( 100.0, 0.0 );
 * // returns 1.0
 *
 * @example
 * var v = pow( 3.141592653589793, 5.0 );
 * // returns ~306.0197
 *
 * @example
 * var v = pow( 3.141592653589793, -0.2 );
 * // returns ~0.7954
 *
 * @example
 * var v = pow( NaN, 3.0 );
 * // returns NaN
 *
 * @example
 * var v = pow( 5.0, NaN );
 * // returns NaN
 *
 * @example
 * var v = pow( NaN, NaN );
 * // returns NaN
 */
 function pow( x, y ) {
 	var ahx; // absolute value high word `x`
 	var ahy; // absolute value high word `y`
 	var ax;  // absolute value `x`
 	var hx;  // high word `x`
 	var lx;  // low word `x`
 	var hy;  // high word `y`
 	var ly;  // low word `y`
 	var sx;  // sign `x`
 	var sy;  // sign `y`
 	var y1;
 	var hp;
 	var lp;
 	var t;
 	var z;   // y prime
 	var j;
 	var i;
 	if ( isnan( x ) || isnan( y ) ) {
 		return NaN;
 	}
 	// Split `y` into high and low words:
 	toWords( WORDS, y );
 	hy = WORDS[ 0 ];
 	ly = WORDS[ 1 ];
 
 	// Special cases `y`...
 	if ( ly === 0 ) {
 		if ( y === 0.0 ) {
 			return 1.0;
 		}
 		if ( y === 1.0 ) {
 			return x;
 		}
 		if ( y === -1.0 ) {
 			return 1.0 / x;
 		}
 		if ( y === 0.5 ) {
 			return sqrt( x );
 		}
 		if ( y === -0.5 ) {
 			return 1.0 / sqrt( x );
 		}
 		if ( y === 2.0 ) {
 			return x * x;
 		}
 		if ( y === 3.0 ) {
 			return x * x * x;
 		}
 		if ( y === 4.0 ) {
 			x *= x;
 			return x * x;
 		}
 		if ( isInfinite( y ) ) {
 			return yIsInfinite( x, y );
 		}
 	}
 	// Split `x` into high and low words:
 	toWords( WORDS, x );
 	hx = WORDS[ 0 ];
 	lx = WORDS[ 1 ];
 
 	// Special cases `x`...
 	if ( lx === 0 ) {
 		if ( hx === 0 ) {
 			return xIsZero( x, y );
 		}
 		if ( x === 1.0 ) {
 			return 1.0;
 		}
 		if (
 			x === -1.0 &&
 			isOdd( y )
 		) {
 			return -1.0;
 		}
 		if ( isInfinite( x ) ) {
 			if ( x === NINF ) {
 				// `pow( 1/x, -y )`
 				return pow( -0.0, -y );
 			}
 			if ( y < 0.0 ) {
 				return 0.0;
 			}
 			return PINF;
 		}
 	}
 	if (
 		x < 0.0 &&
 		isInteger( y ) === false
 	) {
 		// Signal NaN...
 		return (x-x)/(x-x);
 	}
 	ax = abs( x );
 
 	// Remove the sign bits (i.e., get absolute values):
 	ahx = (hx & ABS_MASK)|0; // asm type annotation
 	ahy = (hy & ABS_MASK)|0; // asm type annotation
 
 	// Extract the sign bits:
 	sx = (hx >>> HIGH_NUM_NONSIGN_BITS)|0; // asm type annotation
 	sy = (hy >>> HIGH_NUM_NONSIGN_BITS)|0; // asm type annotation
 
 	// Determine the sign of the result...
 	if ( sx && isOdd( y ) ) {
 		sx = -1.0;
 	} else {
 		sx = 1.0;
 	}
 	// Case 1: `|y|` is huge...
 
 	// |y| > 2^31
 	if ( ahy > HIGH_BIASED_EXP_31 ) {
 		// `|y| > 2^64`, then must over- or underflow...
 		if ( ahy > HIGH_BIASED_EXP_64 ) {
 			return yIsHuge( x, y );
 		}
 		// Over- or underflow if `x` is not close to unity...
 
 		if ( ahx < HIGH_MAX_NEAR_UNITY ) {
 			// y < 0
 			if ( sy === 1 ) {
 				// Signal overflow...
 				return sx * HUGE * HUGE;
 			}
 			// Signal underflow...
 			return sx * TINY * TINY;
 		}
 		if ( ahx > HIGH_BIASED_EXP_0 ) {
 			// y > 0
 			if ( sy === 0 ) {
 				// Signal overflow...
 				return sx * HUGE * HUGE;
 			}
 			// Signal underflow...
 			return sx * TINY * TINY;
 		}
 		// At this point, `|1-x|` is tiny (`<= 2^-20`). Suffice to compute `log(x)` by `x - x^2/2 + x^3/3 - x^4/4`.
 		t = logx( LOG_WORKSPACE, ax );
 	}
 	// Case 2: `|y|` is not huge...
 	else {
 		t = log2ax( LOG_WORKSPACE, ax, ahx );
 	}
 	// Split `y` into `y1 + y2` and compute `(y1+y2) * (t1+t2)`...
 	y1 = setLowWord( y, 0 );
 	lp = ( (y-y1)*t[0] ) + ( y*t[1] );
 	hp = y1 * t[0];
 	z = lp + hp;
 
 	// Note: *can* be more performant to use `getHighWord` and `getLowWord` directly, but using `toWords` looks cleaner.
 	toWords( WORDS, z );
 	j = uint32ToInt32( WORDS[0] );
 	i = uint32ToInt32( WORDS[1] );
 
 	// z >= 1024
 	if ( j >= HIGH_BIASED_EXP_10 ) {
 		// z > 1024
 		if ( ((j-HIGH_BIASED_EXP_10)|i) !== 0 ) {
 			// Signal overflow...
 			return sx * HUGE * HUGE;
 		}
 		if ( (lp+OVT) > (z-hp) ) {
 			// Signal overflow...
 			return sx * HUGE * HUGE;
 		}
 	}
 	// z <= -1075
 	else if ( (j&ABS_MASK) >= HIGH_1075 ) {
 		// z < -1075
 		if ( ((j-HIGH_NEG_1075)|i) !== 0 ) {
 			// signal underflow...
 			return sx * TINY * TINY;
 		}
 		if ( lp <= (z-hp) ) {
 			// signal underflow...
 			return sx * TINY * TINY;
 		}
 	}
 	// Compute `2^(hp+lp)`...
 	z = pow2( j, hp, lp );
 
 	return sx * z;
 }
 
 
 // EXPORTS //
 
 module.exports = pow;