time-to-botec

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

rempio2_medium.js (3276B)

      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/k_rem_pio2.c}. The implementation follows the original, but has been modified for JavaScript.
     22 *
     23 * ```text
     24 * Copyright (C) 1993 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 round = require( './../../../../base/special/round' );
     38 var getHighWord = require( '@stdlib/number/float64/base/get-high-word' );
     39 
     40 
     41 // VARIABLES //
     42 
     43 // 53 bits of 2/π:
     44 var INVPIO2 = 6.36619772367581382433e-01; // 0x3FE45F30, 0x6DC9C883
     45 
     46 // First 33 bits of π/2:
     47 var PIO2_1 = 1.57079632673412561417e+00;  // 0x3FF921FB, 0x54400000
     48 
     49 // PIO2_1T = π/2 - PIO2_1:
     50 var PIO2_1T = 6.07710050650619224932e-11; // 0x3DD0B461, 0x1A626331
     51 
     52 // Another 33 bits of π/2:
     53 var PIO2_2 = 6.07710050630396597660e-11;  // 0x3DD0B461, 0x1A600000
     54 
     55 // PIO2_2T = π/2 - ( PIO2_1 + PIO2_2 ):
     56 var PIO2_2T = 2.02226624879595063154e-21; // 0x3BA3198A, 0x2E037073
     57 
     58 // Another 33 bits of π/2:
     59 var PIO2_3 = 2.02226624871116645580e-21;  // 0x3BA3198A, 0x2E000000
     60 
     61 // PIO2_3T = π/2 - ( PIO2_1 + PIO2_2 + PIO2_3 ):
     62 var PIO2_3T = 8.47842766036889956997e-32; // 0x397B839A, 0x252049C1
     63 
     64 // Exponent mask (2047 => 0x7ff):
     65 var EXPONENT_MASK = 0x7ff|0; // asm type annotation
     66 
     67 
     68 // MAIN //
     69 
     70 /**
     71 * Computes `x - nπ/2 = r` for medium-sized inputs.
     72 *
     73 * @private
     74 * @param {number} x - input value
     75 * @param {uint32} ix - high word of `x`
     76 * @param {(Array|TypedArray|Object)} y - remainder elements
     77 * @returns {integer} factor of `π/2`
     78 */
     79 function rempio2Medium( x, ix, y ) {
     80 	var high;
     81 	var n;
     82 	var t;
     83 	var r;
     84 	var w;
     85 	var i;
     86 	var j;
     87 
     88 	n = round( x * INVPIO2 );
     89 	r = x - ( n * PIO2_1 );
     90 	w = n * PIO2_1T;
     91 
     92 	// First rounding (good to 85 bits)...
     93 	j = (ix >> 20)|0; // asm type annotation
     94 	y[ 0 ] = r - w;
     95 	high = getHighWord( y[0] );
     96 	i = j - ( (high >> 20) & EXPONENT_MASK );
     97 
     98 	// Check if a second iteration is needed (good to 118 bits)...
     99 	if ( i > 16 ) {
    100 		t = r;
    101 		w = n * PIO2_2;
    102 		r = t - w;
    103 		w = (n * PIO2_2T) - ((t-r) - w);
    104 		y[ 0 ] = r - w;
    105 		high = getHighWord( y[0] );
    106 		i = j - ( (high >> 20) & EXPONENT_MASK );
    107 
    108 		// Check if a third iteration is needed (151 bits accumulated)...
    109 		if ( i > 49 ) {
    110 			t = r;
    111 			w = n * PIO2_3;
    112 			r = t - w;
    113 			w = (n * PIO2_3T) - ((t-r) - w);
    114 			y[ 0 ] = r - w;
    115 		}
    116 	}
    117 	y[ 1 ] = (r - y[0]) - w;
    118 	return n;
    119 }
    120 
    121 
    122 // EXPORTS //
    123 
    124 module.exports = rempio2Medium;