kernel_tan.js (5124B)
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, license, and long comment 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_tan.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 setLowWord = require( '@stdlib/number/float64/base/set-low-word' ); 39 var polyvalOdd = require( './polyval_t_odd.js' ); 40 var polyvalEven = require( './polyval_t_even.js' ); 41 42 43 // VARIABLES // 44 45 var PIO4 = 7.85398163397448278999e-01; 46 var PIO4LO = 3.06161699786838301793e-17; 47 var T0 = 3.33333333333334091986e-01; // 3FD55555, 55555563 48 49 // Absolute value mask: 2147483647 => 0x7fffffff => 01111111111111111111111111111111 50 var HIGH_WORD_ABS_MASK = 0x7fffffff|0; // asm type annotation 51 52 53 // MAIN // 54 55 /** 56 * Computes the tangent on \\( \approx\[-\pi/4, \pi/4] \\) (except on -0), \\( \pi/4 \approx 0.7854 \\). 57 * 58 * ## Method 59 * 60 * - Since \\( \tan(-x) = -\tan(x) \\), we need only to consider positive \\( x \\). 61 * 62 * - Callers must return \\( \tan(-0) = -0 \\) without calling here since our odd polynomial is not evaluated in a way that preserves \\( -0 \\). Callers may do the optimization \\( \tan(x) \approx x \\) for tiny \\( x \\). 63 * 64 * - \\( \tan(x) \\) is approximated by a odd polynomial of degree 27 on \\( \[0, 0.67434] \\) 65 * 66 * ```tex 67 * \tan(x) \approx x + T_1 x^3 + \ldots + T_{13} x^{27} 68 * ``` 69 * where 70 * 71 * ```tex 72 * \left| \frac{\tan(x)}{x} - \left( 1 + T_1 x^2 + T_2 x^4 + \ldots + T_{13} x^{26} \right) \right| \le 2^{-59.2} 73 * ``` 74 * 75 * - Note: \\( \tan(x+y) = \tan(x) + \tan'(x) \cdot y \approx \tan(x) + ( 1 + x \cdot x ) \cdot y \\). Therefore, for better accuracy in computing \\( \tan(x+y) \\), let 76 * 77 * ```tex 78 * r = x^3 \cdot \left( T_2+x^2 \cdot (T_3+x^2 \cdot (\ldots+x^2 \cdot (T_{12}+x^2 \cdot T_{13}))) \right) 79 * ``` 80 * 81 * then 82 * 83 * ```tex 84 * \tan(x+y) = x^3 + \left( T_1 \cdot x^2 + (x \cdot (r+y)+y) \right) 85 * ``` 86 * 87 * - For \\( x \\) in \\( \[0.67434, \pi/4] \\), let \\( y = \pi/4 - x \\), then 88 * 89 * ```tex 90 * \tan(x) = \tan\left(\tfrac{\pi}{4}-y\right) = \frac{1-\tan(y)}{1+\tan(y)} \\ 91 * = 1 - 2 \cdot \left( \tan(y) - \tfrac{\tan(y)^2}{1+\tan(y)} \right) 92 * ``` 93 * 94 * 95 * @param {number} x - input value (in radians, assumed to be bounded by ~π/4 in magnitude) 96 * @param {number} y - tail of `x` 97 * @param {integer} k - indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned 98 * @returns {number} tangent 99 * 100 * @example 101 * var out = kernelTan( 3.141592653589793/4.0, 0.0, 1 ); 102 * // returns ~1.0 103 * 104 * @example 105 * var out = kernelTan( 3.141592653589793/4.0, 0.0, -1 ); 106 * // returns ~-1.0 107 * 108 * @example 109 * var out = kernelTan( 3.141592653589793/6.0, 0.0, 1 ); 110 * // returns ~0.577 111 * 112 * @example 113 * var out = kernelTan( 0.664, 5.288e-17, 1 ); 114 * // returns ~0.783 115 * 116 * @example 117 * var out = kernelTan( NaN, 0.0, 1 ); 118 * // returns NaN 119 * 120 * @example 121 * var out = kernelTan( 3.0, NaN, 1 ); 122 * // returns NaN 123 * 124 * @example 125 * var out = kernelTan( NaN, NaN, 1 ); 126 * // returns NaN 127 */ 128 function kernelTan( x, y, k ) { 129 var hx; 130 var ix; 131 var a; 132 var r; 133 var s; 134 var t; 135 var v; 136 var w; 137 var z; 138 139 hx = getHighWord( x ); 140 141 // High word of |x|: 142 ix = (hx & HIGH_WORD_ABS_MASK)|0; // asm type annotation 143 144 // Case: |x| >= 0.6744 145 if ( ix >= 0x3FE59428 ) { 146 if ( x < 0 ) { 147 x = -x; 148 y = -y; 149 } 150 z = PIO4 - x; 151 w = PIO4LO - y; 152 x = z + w; 153 y = 0.0; 154 } 155 z = x * x; 156 w = z * z; 157 158 // Break x^5*(T[1]+x^2*T[2]+...) into x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + x^5(x^2*(T[2]+x^4*T[4]+...+x^22*T[12]))... 159 r = polyvalOdd( w ); 160 v = z * polyvalEven( w ); 161 s = z * x; 162 r = y + (z * ((s * (r + v)) + y)); 163 r += T0 * s; 164 w = x + r; 165 if ( ix >= 0x3FE59428 ) { 166 v = k; 167 return ( 1.0 - ( (hx >> 30) & 2 ) ) * ( v - (2.0 * (x - ((w * w / (w + v)) - r)) )); // eslint-disable-line max-len 168 } 169 if ( k === 1 ) { 170 return w; 171 } 172 // Compute -1/(x+r) accurately... 173 z = w; 174 setLowWord( z, 0 ); 175 v = r - (z - x); // z + v = r + x 176 a = -1.0 / w; // a = -1/w 177 t = a; 178 setLowWord( t, 0 ); 179 s = 1.0 + (t * z); 180 return t + (a * (s + (t * v))); 181 } 182 183 184 // EXPORTS // 185 186 module.exports = kernelTan;