sici.js (4918B)
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 original C code, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. The implementation follows the original, but has been modified for JavaScript. 22 * 23 * ```text 24 * Copyright 1984, 1987, 1989 by Stephen L. Moshier 25 * 26 * Some software in this archive may be from the book _Methods and Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster International, 1989) or from the Cephes Mathematical Library, a commercial product. In either event, it is copyrighted by the author. What you see here may be used freely but it comes with no support or guarantee. 27 * 28 * Stephen L. Moshier 29 * moshier@na-net.ornl.gov 30 * ``` 31 */ 32 33 'use strict'; 34 35 // MODULES // 36 37 var isInfinite = require( './../../../../base/assert/is-infinite' ); 38 var isnan = require( './../../../../base/assert/is-nan' ); 39 var cos = require( './../../../../base/special/cos' ); 40 var sin = require( './../../../../base/special/sin' ); 41 var ln = require( './../../../../base/special/ln' ); 42 var HALF_PI = require( '@stdlib/constants/float64/half-pi' ); 43 var GAMMA = require( '@stdlib/constants/float64/eulergamma' ); 44 var NINF = require( '@stdlib/constants/float64/ninf' ); 45 var polyvalFN4 = require( './polyval_fn4.js' ); 46 var polyvalFD4 = require( './polyval_fd4.js' ); 47 var polyvalFN8 = require( './polyval_fn8.js' ); 48 var polyvalFD8 = require( './polyval_fd8.js' ); 49 var polyvalGN4 = require( './polyval_gn4.js' ); 50 var polyvalGD4 = require( './polyval_gd4.js' ); 51 var polyvalGN8 = require( './polyval_gn8.js' ); 52 var polyvalGD8 = require( './polyval_gd8.js' ); 53 var polyvalSN = require( './polyval_sn.js' ); 54 var polyvalSD = require( './polyval_sd.js' ); 55 var polyvalCN = require( './polyval_cn.js' ); 56 var polyvalCD = require( './polyval_cd.js' ); 57 58 59 // MAIN // 60 61 /** 62 * Computes the sine and cosine integrals. 63 * 64 * ## Method 65 * 66 * - The integrals are approximated by rational functions. 67 * 68 * - For \\( x > 8 \\), auxiliary functions \\( f(x) \\) and \\( g(x) \\) are employed such that 69 * 70 * ```tex 71 * \operatorname{Ci}(x) = f(x) \sin(x) - g(x) \cos(x) \\ 72 * \operatorname{Si}(x) = \pi/2 - f(x) \cos(x) - g(x) \sin(x) 73 * ``` 74 * 75 * ## Notes 76 * 77 * - Absolute error on test interval \\( \[0,50\] \\), except relative when greater than \\( 1 \\): 78 * 79 * | arithmetic | function | # trials | peak | rms | 80 * |:----------:|:-----------:|:--------:|:-------:|:-------:| 81 * | IEEE | Si | 30000 | 4.4e-16 | 7.3e-17 | 82 * | IEEE | Ci | 30000 | 6.9e-16 | 5.1e-17 | 83 * 84 * 85 * @private 86 * @param {(Array|TypedArray|Object)} out - output array 87 * @param {number} x - input value 88 * @returns {(Array|TypedArray|Object)} output array 89 * 90 * @example 91 * var v = sici( [ 0.0, 0.0 ], 3.0 ); 92 * // returns [ ~1.849, ~0.12 ] 93 * 94 * @example 95 * var v = sici( [ 0.0, 0.0 ], 0.0 ); 96 * // returns [ 0.0, -Infinity ] 97 * 98 * @example 99 * var v = sici( [ 0.0, 0.0 ], -9.0 ); 100 * // returns [ ~-1.665, ~0.055 ] 101 * 102 * @example 103 * var v = sici( [ 0.0, 0.0 ], NaN ); 104 * // returns [ NaN, NaN ] 105 */ 106 function sici( out, x ) { 107 var sgn; 108 var si; 109 var ci; 110 var c; 111 var f; 112 var g; 113 var s; 114 var z; 115 116 if ( isnan( x ) ) { 117 out[ 0 ] = NaN; 118 out[ 1 ] = NaN; 119 return out; 120 } 121 if ( x < 0.0 ) { 122 sgn = -1; 123 x = -x; 124 } else { 125 sgn = 0; 126 } 127 if ( x === 0.0 ) { 128 out[ 0 ] = 0.0; 129 out[ 1 ] = NINF; 130 return out; 131 } 132 if ( x > 1.0e9 ) { 133 if ( isInfinite( x ) ) { 134 if ( sgn === -1 ) { 135 si = -HALF_PI; 136 ci = NaN; 137 } else { 138 si = HALF_PI; 139 ci = 0.0; 140 } 141 out[ 0 ] = si; 142 out[ 1 ] = ci; 143 return out; 144 } 145 si = HALF_PI - ( cos( x ) / x ); 146 ci = sin( x ) / x; 147 } 148 if ( x > 4.0 ) { 149 s = sin( x ); 150 c = cos( x ); 151 z = 1.0 / ( x*x ); 152 if ( x < 8.0 ) { 153 f = polyvalFN4( z ) / ( x * polyvalFD4( z ) ); 154 g = z * polyvalGN4( z ) / polyvalGD4( z ); 155 } else { 156 f = polyvalFN8( z ) / ( x * polyvalFD8( z ) ); 157 g = z * polyvalGN8( z ) / polyvalGD8( z ); 158 } 159 si = HALF_PI - ( f*c ) - ( g*s ); 160 if ( sgn ) { 161 si = -si; 162 } 163 ci = ( f*s ) - ( g*c ); 164 out[ 0 ] = si; 165 out[ 1 ] = ci; 166 return out; 167 } 168 z = x * x; 169 s = x * polyvalSN( z ) / polyvalSD( z ); 170 c = z * polyvalCN( z ) / polyvalCD( z ); 171 if ( sgn ) { 172 s = -s; 173 } 174 si = s; 175 ci = GAMMA + ln( x ) + c; // real part if x < 0 176 out[ 0 ] = si; 177 out[ 1 ] = ci; 178 return out; 179 } 180 181 182 // EXPORTS // 183 184 module.exports = sici;