rempio2.js (7275B)
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/e_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 * Optimized by Bruce D. Evans. 32 * ``` 33 */ 34 35 'use strict'; 36 37 // MODULES // 38 39 var getHighWord = require( '@stdlib/number/float64/base/get-high-word' ); 40 var getLowWord = require( '@stdlib/number/float64/base/get-low-word' ); 41 var fromWords = require( '@stdlib/number/float64/base/from-words' ); 42 var rempio2Kernel = require( './kernel_rempio2.js' ); 43 var rempio2Medium = require( './rempio2_medium.js' ); 44 45 46 // VARIABLES // 47 48 var ZERO = 0.00000000000000000000e+00; // 0x00000000, 0x00000000 49 var TWO24 = 1.67772160000000000000e+07; // 0x41700000, 0x00000000 50 51 // 33 bits of π/2: 52 var PIO2_1 = 1.57079632673412561417e+00; // 0x3FF921FB, 0x54400000 53 54 // PIO2_1T = π/2 - PIO2_1: 55 var PIO2_1T = 6.07710050650619224932e-11; // 0x3DD0B461, 0x1A626331 56 var TWO_PIO2_1T = 2.0 * PIO2_1T; 57 var THREE_PIO2_1T = 3.0 * PIO2_1T; 58 var FOUR_PIO2_1T = 4.0 * PIO2_1T; 59 60 // Absolute value mask: 0x7fffffff = 2147483647 => 01111111111111111111111111111111 61 var ABS_MASK = 0x7fffffff|0; // asm type annotation 62 63 // Exponent mask: 0x7ff00000 = 2146435072 => 01111111111100000000000000000000 64 var EXPONENT_MASK = 0x7ff00000|0; // asm type annotation 65 66 // High word significand mask: 0xfffff = 1048575 => 00000000000011111111111111111111 67 var SIGNIFICAND_MASK = 0xfffff|0; // asm type annotation 68 69 // High word significand for π and π/2: 0x921fb = 598523 => 00000000000010010010000111111011 70 var PI_HIGH_WORD_SIGNIFICAND = 0x921fb|0; // asm type annotation 71 72 // High word for π/4: 0x3fe921fb = 1072243195 => 00111111111010010010000111111011 73 var PIO4_HIGH_WORD = 0x3fe921fb|0; // asm type annotation 74 75 // High word for 3π/4: 0x4002d97c = 1073928572 => 01000000000000101101100101111100 76 var THREE_PIO4_HIGH_WORD = 0x4002d97c|0; // asm type annotation 77 78 // High word for 5π/4: 0x400f6a7a = 1074752122 => 01000000000011110110101001111010 79 var FIVE_PIO4_HIGH_WORD = 0x400f6a7a|0; // asm type annotation 80 81 // High word for 6π/4: 0x4012d97c = 1074977148 => 01000000000100101101100101111100 82 var THREE_PIO2_HIGH_WORD = 0x4012d97c|0; // asm type annotation 83 84 // High word for 7π/4: 0x4015fdbc = 1075183036 => 01000000000101011111110110111100 85 var SEVEN_PIO4_HIGH_WORD = 0x4015fdbc|0; // asm type annotation 86 87 // High word for 8π/4: 0x401921fb = 1075388923 => 01000000000110010010000111111011 88 var TWO_PI_HIGH_WORD = 0x401921fb|0; // asm type annotation 89 90 // High word for 9π/4: 0x401c463b = 1075594811 => 01000000000111000100011000111011 91 var NINE_PIO4_HIGH_WORD = 0x401c463b|0; // asm type annotation 92 93 // 2^20*π/2 = 1647099.3291652855 => 0100000100111001001000011111101101010100010001000010110100011000 => high word => 0x413921fb = 1094263291 => 01000001001110010010000111111011 94 var MEDIUM = 0x413921fb|0; // asm type annotation 95 96 // Arrays for storing temporary values: 97 var TX = [ 0.0, 0.0, 0.0 ]; // WARNING: not thread safe 98 var TY = [ 0.0, 0.0 ]; // WARNING: not thread safe 99 100 101 // MAIN // 102 103 /** 104 * Computes `x - nπ/2 = r`. 105 * 106 * ## Notes 107 * 108 * - Returns `n` and stores the remainder `r` as two numbers `y[0]` and `y[1]`, such that `y[0]+y[1] = r`. 109 * 110 * 111 * @param {number} x - input value 112 * @param {(Array|TypedArray|Object)} y - remainder elements 113 * @returns {integer} factor of `π/2` 114 * 115 * @example 116 * var y = [ 0.0, 0.0 ]; 117 * var n = rempio2( 128.0, y ); 118 * // returns 81 119 * 120 * var y1 = y[ 0 ]; 121 * // returns ~0.765 122 * 123 * var y2 = y[ 1 ]; 124 * // returns ~3.618e-17 125 * 126 * @example 127 * var y = [ 0.0, 0.0 ]; 128 * var n = rempio2( NaN, y ); 129 * // returns 0 130 * 131 * var y1 = y[ 0 ]; 132 * // returns NaN 133 * 134 * var y2 = y[ 1 ]; 135 * // returns NaN 136 */ 137 function rempio2( x, y ) { 138 var low; 139 var e0; 140 var hx; 141 var ix; 142 var nx; 143 var i; 144 var n; 145 var z; 146 147 hx = getHighWord( x ); 148 ix = (hx & ABS_MASK)|0; // asm type annotation 149 150 // Case: |x| ~<= π/4 (no need for reduction) 151 if ( ix <= PIO4_HIGH_WORD ) { 152 y[ 0 ] = x; 153 y[ 1 ] = 0.0; 154 return 0; 155 } 156 // Case: |x| ~<= 5π/4 157 if ( ix <= FIVE_PIO4_HIGH_WORD ) { 158 // Case: |x| ~= π/2 or π 159 if ( (ix & SIGNIFICAND_MASK) === PI_HIGH_WORD_SIGNIFICAND ) { 160 // Cancellation => use medium case 161 return rempio2Medium( x, ix, y ); 162 } 163 // Case: |x| ~<= 3π/4 164 if ( ix <= THREE_PIO4_HIGH_WORD ) { 165 if ( x > 0.0 ) { 166 z = x - PIO2_1; 167 y[ 0 ] = z - PIO2_1T; 168 y[ 1 ] = (z - y[0]) - PIO2_1T; 169 return 1; 170 } 171 z = x + PIO2_1; 172 y[ 0 ] = z + PIO2_1T; 173 y[ 1 ] = (z - y[0]) + PIO2_1T; 174 return -1; 175 } 176 if ( x > 0.0 ) { 177 z = x - ( 2.0*PIO2_1 ); 178 y[ 0 ] = z - TWO_PIO2_1T; 179 y[ 1 ] = (z - y[0]) - TWO_PIO2_1T; 180 return 2; 181 } 182 z = x + ( 2.0*PIO2_1 ); 183 y[ 0 ] = z + TWO_PIO2_1T; 184 y[ 1 ] = (z - y[0]) + TWO_PIO2_1T; 185 return -2; 186 } 187 // Case: |x| ~<= 9π/4 188 if ( ix <= NINE_PIO4_HIGH_WORD ) { 189 // Case: |x| ~<= 7π/4 190 if ( ix <= SEVEN_PIO4_HIGH_WORD ) { 191 // Case: |x| ~= 3π/2 192 if ( ix === THREE_PIO2_HIGH_WORD ) { 193 return rempio2Medium( x, ix, y ); 194 } 195 if ( x > 0.0 ) { 196 z = x - ( 3.0*PIO2_1 ); 197 y[ 0 ] = z - THREE_PIO2_1T; 198 y[ 1 ] = (z - y[0]) - THREE_PIO2_1T; 199 return 3; 200 } 201 z = x + ( 3.0*PIO2_1 ); 202 y[ 0 ] = z + THREE_PIO2_1T; 203 y[ 1 ] = (z - y[0]) + THREE_PIO2_1T; 204 return -3; 205 } 206 // Case: |x| ~= 4π/2 207 if ( ix === TWO_PI_HIGH_WORD ) { 208 return rempio2Medium( x, ix, y ); 209 } 210 if ( x > 0.0 ) { 211 z = x - ( 4.0*PIO2_1 ); 212 y[ 0 ] = z - FOUR_PIO2_1T; 213 y[ 1 ] = (z - y[0]) - FOUR_PIO2_1T; 214 return 4; 215 } 216 z = x + ( 4.0*PIO2_1 ); 217 y[ 0 ] = z + FOUR_PIO2_1T; 218 y[ 1 ] = (z - y[0]) + FOUR_PIO2_1T; 219 return -4; 220 } 221 // Case: |x| ~< 2^20*π/2 (medium size) 222 if ( ix < MEDIUM ) { 223 return rempio2Medium( x, ix, y ); 224 } 225 // Case: x is NaN or infinity 226 if ( ix >= EXPONENT_MASK ) { 227 y[ 0 ] = NaN; 228 y[ 1 ] = NaN; 229 return 0.0; 230 } 231 // Set z = scalbn(|x|, ilogb(x)-23)... 232 low = getLowWord( x ); 233 e0 = (ix >> 20) - 1046; // `e0 = ilogb(z) - 23` => unbiased exponent minus 23 234 z = fromWords( ix - ((e0 << 20)|0), low ); 235 for ( i = 0; i < 2; i++ ) { 236 TX[ i ] = z|0; 237 z = (z - TX[i]) * TWO24; 238 } 239 TX[ 2 ] = z; 240 nx = 3; 241 while ( TX[ nx-1 ] === ZERO ) { 242 // Skip zero term... 243 nx -= 1; 244 } 245 n = rempio2Kernel( TX, TY, e0, nx, 1 ); 246 if ( x < 0.0 ) { 247 y[ 0 ] = -TY[ 0 ]; 248 y[ 1 ] = -TY[ 1 ]; 249 return -n; 250 } 251 y[ 0 ] = TY[ 0 ]; 252 y[ 1 ] = TY[ 1 ]; 253 return n; 254 } 255 256 257 // EXPORTS // 258 259 module.exports = rempio2;