time-to-botec

Benchmark sampling in different programming languages
Log | Files | Refs | README

pow2.js (4246B)

      1 /**
      2 * @license Apache-2.0
      3 *
      4 * Copyright (c) 2018 The Stdlib Authors.
      5 *
      6 * Licensed under the Apache License, Version 2.0 (the "License");
      7 * you may not use this file except in compliance with the License.
      8 * You may obtain a copy of the License at
      9 *
     10 *    http://www.apache.org/licenses/LICENSE-2.0
     11 *
     12 * Unless required by applicable law or agreed to in writing, software
     13 * distributed under the License is distributed on an "AS IS" BASIS,
     14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     15 * See the License for the specific language governing permissions and
     16 * limitations under the License.
     17 *
     18 *
     19 * ## Notice
     20 *
     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/s_pow.c}. The implementation follows the original, but has been modified for JavaScript.
     22 *
     23 * ```text
     24 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
     25 *
     26 * Developed at SunPro, a Sun Microsystems, Inc. business.
     27 * Permission to use, copy, modify, and distribute this
     28 * software is freely granted, provided that this notice
     29 * is preserved.
     30 * ```
     31 */
     32 
     33 'use strict';
     34 
     35 // MODULES //
     36 
     37 var getHighWord = require( '@stdlib/number/float64/base/get-high-word' );
     38 var setHighWord = require( '@stdlib/number/float64/base/set-high-word' );
     39 var setLowWord = require( '@stdlib/number/float64/base/set-low-word' );
     40 var uint32ToInt32 = require( '@stdlib/number/uint32/base/to-int32' );
     41 var ldexp = require( './../../../../base/special/ldexp' );
     42 var LN2 = require( '@stdlib/constants/float64/ln-two' );
     43 var BIAS = require( '@stdlib/constants/float64/exponent-bias' );
     44 var polyvalP = require( './polyval_p.js' );
     45 
     46 
     47 // VARIABLES //
     48 
     49 // 0x7fffffff = 2147483647 => 0 11111111111 11111111111111111111
     50 var ABS_MASK = 0x7fffffff|0; // asm type annotation
     51 
     52 // 0x000fffff = 1048575 => 0 00000000000 11111111111111111111
     53 var HIGH_SIGNIFICAND_MASK = 0x000fffff|0; // asm type annotation
     54 
     55 // 0x00100000 = 1048576 => 0 00000000001 00000000000000000000 => biased exponent: 1 = -1022+1023 => 2^-1022
     56 var HIGH_MIN_NORMAL_EXP = 0x00100000|0; // asm type annotation
     57 
     58 // 0x3fe00000 = 1071644672 => 0 01111111110 00000000000000000000 => biased exponent: 1022 = -1+1023 => 2^-1
     59 var HIGH_BIASED_EXP_NEG_1 = 0x3fe00000|0; // asm type annotation
     60 
     61 // TODO: consider making into an external constant
     62 var HIGH_NUM_SIGNIFICAND_BITS = 20|0; // asm type annotation
     63 
     64 // High: LN2
     65 var LN2_HI = 6.93147182464599609375e-01; // 0x3FE62E43, 0x00000000
     66 
     67 // Low: LN2
     68 var LN2_LO = -1.90465429995776804525e-09; // 0xBE205C61, 0x0CA86C39
     69 
     70 
     71 // MAIN //
     72 
     73 /**
     74 * Computes \\(2^{\mathrm{hp} + \mathrm{lp}\\).
     75 *
     76 * @private
     77 * @param {number} j - high word of `hp + lp`
     78 * @param {number} hp - first power summand
     79 * @param {number} lp - second power summand
     80 * @returns {number} function value
     81 *
     82 * @example
     83 * var z = pow2( 1065961648, -0.3398475646972656, -0.000002438187359100815 );
     84 * // returns ~0.79
     85 */
     86 function pow2( j, hp, lp ) {
     87 	var tmp;
     88 	var t1;
     89 	var t;
     90 	var r;
     91 	var u;
     92 	var v;
     93 	var w;
     94 	var z;
     95 	var n;
     96 	var i;
     97 	var k;
     98 
     99 	i = (j & ABS_MASK)|0; // asm type annotation
    100 	k = ((i>>HIGH_NUM_SIGNIFICAND_BITS) - BIAS)|0; // asm type annotation
    101 	n = 0;
    102 
    103 	// `|z| > 0.5`, set `n = z+0.5`
    104 	if ( i > HIGH_BIASED_EXP_NEG_1 ) {
    105 		n = (j + (HIGH_MIN_NORMAL_EXP>>(k+1)))>>>0; // asm type annotation
    106 		k = (((n & ABS_MASK)>>HIGH_NUM_SIGNIFICAND_BITS) - BIAS)|0; // new k for n
    107 		tmp = ((n & ~(HIGH_SIGNIFICAND_MASK >> k)))>>>0; // asm type annotation
    108 		t = setHighWord( 0.0, tmp );
    109 		n = (((n & HIGH_SIGNIFICAND_MASK)|HIGH_MIN_NORMAL_EXP) >> (HIGH_NUM_SIGNIFICAND_BITS-k))>>>0; // eslint-disable-line max-len
    110 		if ( j < 0 ) {
    111 			n = -n;
    112 		}
    113 		hp -= t;
    114 	}
    115 	t = lp + hp;
    116 	t = setLowWord( t, 0 );
    117 	u = t * LN2_HI;
    118 	v = ( (lp - (t-hp))*LN2 ) + ( t*LN2_LO );
    119 	z = u + v;
    120 	w = v - (z - u);
    121 	t = z * z;
    122 	t1 = z - ( t*polyvalP( t ) );
    123 	r = ( (z*t1) / (t1-2.0) ) - ( w + (z*w) );
    124 	z = 1.0 - (r - z);
    125 	j = getHighWord( z );
    126 	j = uint32ToInt32( j );
    127 	j += (n << HIGH_NUM_SIGNIFICAND_BITS)>>>0; // asm type annotation
    128 
    129 	// Check for subnormal output...
    130 	if ( (j>>HIGH_NUM_SIGNIFICAND_BITS) <= 0 ) {
    131 		z = ldexp( z, n );
    132 	} else {
    133 		z = setHighWord( z, j );
    134 	}
    135 	return z;
    136 }
    137 
    138 
    139 // EXPORTS //
    140 
    141 module.exports = pow2;