time-to-botec

kernel_rempio2.js (8937B)

---

 /**
 * @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/k_rem_pio2.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.
 * ```
 */
 
 /* eslint-disable array-element-newline */
 
 'use strict';
 
 // MODULES //
 
 var floor = require( './../../../../base/special/floor' );
 var ldexp = require( './../../../../base/special/ldexp' );
 
 
 // VARIABLES //
 
 /*
 * Table of constants for `2/π` (`396` hex digits, `476` decimal).
 *
 * Integer array which contains the (`24*i`)-th to (`24*i+23`)-th bit of `2/π` after binary point. The corresponding floating value is
 *
 * ```tex
 * \operatorname{ipio2}[i] \cdot 2^{-24(i+1)}
 * ```
 *
 * This table must have at least `(e0-3)/24 + jk` terms. For quad precision (`e0 <= 16360`, `jk = 6`), this is `686`.
 */
 var IPIO2 = [
 	0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
 	0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
 	0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
 	0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
 	0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
 	0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
 	0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
 	0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
 	0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
 	0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
 	0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B
 ];
 
 // Double precision array, obtained by cutting `π/2` into `24` bits chunks...
 var PIO2 = [
 	1.57079625129699707031e+00, // 0x3FF921FB, 0x40000000
 	7.54978941586159635335e-08, // 0x3E74442D, 0x00000000
 	5.39030252995776476554e-15, // 0x3CF84698, 0x80000000
 	3.28200341580791294123e-22, // 0x3B78CC51, 0x60000000
 	1.27065575308067607349e-29, // 0x39F01B83, 0x80000000
 	1.22933308981111328932e-36, // 0x387A2520, 0x40000000
 	2.73370053816464559624e-44, // 0x36E38222, 0x80000000
 	2.16741683877804819444e-51  // 0x3569F31D, 0x00000000
 ];
 var TWO24 = 1.67772160000000000000e+07;  // 0x41700000, 0x00000000
 var TWON24 = 5.96046447753906250000e-08; // 0x3E700000, 0x00000000
 
 // Arrays for storing temporary values (note that, in C, this is not thread safe):
 var F = zeros( 20 );
 var Q = zeros( 20 );
 var FQ = zeros( 20 );
 var IQ = zeros( 20 );
 
 
 // FUNCTIONS //
 
 /**
 * Returns an array of zeros.
 *
 * @private
 * @param {NonNegativeInteger} len - array length
 * @returns {NonNegativeIntegerArray} output array
 */
 function zeros( len ) {
 	var out;
 	var i;
 
 	out = [];
 	for ( i = 0; i < len; i++ ) {
 		out.push( 0.0 );
 	}
 	return out;
 }
 
 /**
 * Performs the computation for `kernelRempio2()`.
 *
 * @private
 * @param {PositiveNumber} x - input value
 * @param {(Array|TypedArray|Object)} y - output object for storing double precision numbers
 * @param {integer} jz - number of terms of `ipio2[]` used
 * @param {Array<integer>} q - array with integral values, representing the 24-bits chunk of the product of `x` and `2/π`
 * @param {integer} q0 - the corresponding exponent of `q[0]` (the exponent for `q[i]` would be `q0-24*i`)
 * @param {integer} jk - `jk+1` is the initial number of terms of `IPIO2[]` needed in the computation
 * @param {integer} jv - index for pointing to the suitable `ipio2[]` for the computation
 * @param {integer} jx - `nx - 1`
 * @param {Array<number>} f - `IPIO2[]` in floating point
 * @returns {number} last three binary digits of `N`
 */
 function compute( x, y, jz, q, q0, jk, jv, jx, f ) {
 	var carry;
 	var fw;
 	var ih;
 	var jp;
 	var i;
 	var k;
 	var n;
 	var j;
 	var z;
 
 	// `jp+1` is the number of terms in `PIO2[]` needed:
 	jp = jk;
 
 	// Distill `q[]` into `IQ[]` in reverse order...
 	z = q[ jz ];
 	j = jz;
 	for ( i = 0; j > 0; i++ ) {
 		fw = ( TWON24 * z )|0;
 		IQ[ i ] = ( z - (TWO24*fw) )|0;
 		z = q[ j-1 ] + fw;
 		j -= 1;
 	}
 	// Compute `n`...
 	z = ldexp( z, q0 );
 	z -= 8.0 * floor( z*0.125 ); // Trim off integer >= 8
 	n = z|0;
 	z -= n;
 	ih = 0;
 	if ( q0 > 0 ) {
 		// Need `IQ[jz-1]` to determine `n`...
 		i = ( IQ[ jz-1 ] >> (24-q0) );
 		n += i;
 		IQ[ jz-1 ] -= ( i << (24-q0) );
 		ih = ( IQ[ jz-1 ] >> (23-q0) );
 	}
 	else if ( q0 === 0 ) {
 		ih = ( IQ[ jz-1 ] >> 23 );
 	}
 	else if ( z >= 0.5 ) {
 		ih = 2;
 	}
 	// Case: q > 0.5
 	if ( ih > 0 ) {
 		n += 1;
 		carry = 0;
 
 		// Compute `1-q`:
 		for ( i = 0; i < jz; i++ ) {
 			j = IQ[ i ];
 			if ( carry === 0 ) {
 				if ( j !== 0 ) {
 					carry = 1;
 					IQ[ i ] = 0x1000000 - j;
 				}
 			} else {
 				IQ[ i ] = 0xffffff - j;
 			}
 		}
 		if ( q0 > 0 ) {
 			// Rare case: chance is 1 in 12...
 			switch ( q0 ) { // eslint-disable-line default-case
 			case 1:
 				IQ[ jz-1 ] &= 0x7fffff;
 				break;
 			case 2:
 				IQ[ jz-1 ] &= 0x3fffff;
 				break;
 			}
 		}
 		if ( ih === 2 ) {
 			z = 1.0 - z;
 			if ( carry !== 0 ) {
 				z -= ldexp( 1.0, q0 );
 			}
 		}
 	}
 	// Check if re-computation is needed...
 	if ( z === 0.0 ) {
 		j = 0;
 		for ( i = jz-1; i >= jk; i-- ) {
 			j |= IQ[ i ];
 		}
 		if ( j === 0 ) {
 			// Need re-computation...
 			for ( k = 1; IQ[ jk-k ] === 0; k++ ) {
 				// `k` is the number of terms needed...
 			}
 			for ( i = jz+1; i <= jz+k; i++ ) {
 				// Add `q[jz+1]` to `q[jz+k]`...
 				f[ jx+i ] = IPIO2[ jv+i ];
 				fw = 0.0;
 				for ( j = 0; j <= jx; j++ ) {
 					fw += x[ j ] * f[ jx + (i-j) ];
 				}
 				q[ i ] = fw;
 			}
 			jz += k;
 			return compute( x, y, jz, q, q0, jk, jv, jx, f );
 		}
 	}
 	// Chop off zero terms...
 	if ( z === 0.0 ) {
 		jz -= 1;
 		q0 -= 24;
 		while ( IQ[ jz ] === 0 ) {
 			jz -= 1;
 			q0 -= 24;
 		}
 	} else {
 		// Break `z` into 24-bit if necessary...
 		z = ldexp( z, -q0 );
 		if ( z >= TWO24 ) {
 			fw = (TWON24*z)|0;
 			IQ[ jz ] = ( z - (TWO24*fw) )|0;
 			jz += 1;
 			q0 += 24;
 			IQ[ jz ] = fw;
 		} else {
 			IQ[ jz ] = z|0;
 		}
 	}
 	// Convert integer "bit" chunk to floating-point value...
 	fw = ldexp( 1.0, q0 );
 	for ( i = jz; i >= 0; i-- ) {
 		q[ i ] = fw * IQ[i];
 		fw *= TWON24;
 	}
 	// Compute `PIO2[0,...,jp]*q[jz,...,0]`...
 	for ( i = jz; i >= 0; i-- ) {
 		fw = 0.0;
 		for ( k = 0; k <= jp && k <= jz-i; k++ ) {
 			fw += PIO2[ k ] * q[ i+k ];
 		}
 		FQ[ jz-i ] = fw;
 	}
 	// Compress `FQ[]` into `y[]`...
 	fw = 0.0;
 	for ( i = jz; i >= 0; i-- ) {
 		fw += FQ[ i ];
 	}
 	if ( ih === 0 ) {
 		y[ 0 ] = fw;
 	} else {
 		y[ 0 ] = -fw;
 	}
 	fw = FQ[ 0 ] - fw;
 	for ( i = 1; i <= jz; i++ ) {
 		fw += FQ[i];
 	}
 	if ( ih === 0 ) {
 		y[ 1 ] = fw;
 	} else {
 		y[ 1 ] = -fw;
 	}
 	return ( n & 7 );
 }
 
 
 // MAIN //
 
 /**
 * Returns the last three binary digits of `N` with `y = x - Nπ/2` so that `|y| < π/2`.
 *
 * ## Method
 *
 * -   The method is to compute the integer (`mod 8`) and fraction parts of `2x/π` without doing the full multiplication. In general, we skip the part of the product that is known to be a huge integer (more accurately, equals `0 mod 8` ). Thus, the number of operations is independent of the exponent of the input.
 *
 * @private
 * @param {PositiveNumber} x - input value
 * @param {(Array|TypedArray|Object)} y - remainder elements
 * @param {PositiveInteger} e0 - the exponent of `x[0]` (must be <= 16360)
 * @param {PositiveInteger} nx - dimension of `x[]`
 * @returns {number} last three binary digits of `N`
 */
 function kernelRempio2( x, y, e0, nx ) {
 	var fw;
 	var jk;
 	var jv;
 	var jx;
 	var jz;
 	var q0;
 	var i;
 	var j;
 	var m;
 
 	// Initialize `jk` for double-precision floating-point numbers:
 	jk = 4;
 
 	// Determine `jx`, `jv`, `q0` (note that `q0 < 3`):
 	jx = nx - 1;
 	jv = ( (e0 - 3) / 24 )|0;
 	if ( jv < 0 ) {
 		jv = 0;
 	}
 	q0 = e0 - (24 * (jv + 1));
 
 	// Set up `F[0]` to `F[jx+jk]` where `F[jx+jk] = IPIO2[jv+jk]`:
 	j = jv - jx;
 	m = jx + jk;
 	for ( i = 0; i <= m; i++ ) {
 		if ( j < 0 ) {
 			F[ i ] = 0.0;
 		} else {
 			F[ i ] = IPIO2[ j ];
 		}
 		j += 1;
 	}
 	// Compute `Q[0],Q[1],...,Q[jk]`:
 	for ( i = 0; i <= jk; i++ ) {
 		fw = 0.0;
 		for ( j = 0; j <= jx; j++ ) {
 			fw += x[ j ] * F[ jx + (i-j) ];
 		}
 		Q[ i ] = fw;
 	}
 	jz = jk;
 	return compute( x, y, jz, Q, q0, jk, jv, jx, F );
 }
 
 
 // EXPORTS //
 
 module.exports = kernelRempio2;