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