time-to-botec

Benchmark sampling in different programming languages
Log | Files | Refs | README

y1.js (3909B)

      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 and copyright notice are from the [Boost library]{@link https://github.com/boostorg/math/blob/develop/include/boost/math/special_functions/detail/bessel_y1.hpp}. The implementation has been modified for JavaScript.
     22 *
     23 * ```text
     24 * Copyright Xiaogang Zhang, 2006.
     25 *
     26 * Use, modification and distribution are subject to the
     27 * Boost Software License, Version 1.0. (See accompanying file
     28 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
     29 * ```
     30 */
     31 
     32 'use strict';
     33 
     34 // MODULES //
     35 
     36 var ln = require( './../../../../base/special/ln' );
     37 var sqrt = require( './../../../../base/special/sqrt' );
     38 var PI = require( '@stdlib/constants/float64/pi' );
     39 var SQRT_PI = require( '@stdlib/constants/float64/sqrt-pi' );
     40 var NINF = require( '@stdlib/constants/float64/ninf' );
     41 var PINF = require( '@stdlib/constants/float64/pinf' );
     42 var sincos = require( './../../../../base/special/sincos' );
     43 var besselj1 = require( './../../../../base/special/besselj1' );
     44 var poly1 = require( './rational_p1q1.js' );
     45 var poly2 = require( './rational_p2q2.js' );
     46 var polyC = require( './rational_pcqc.js' );
     47 var polyS = require( './rational_psqs.js' );
     48 
     49 
     50 // VARIABLES //
     51 
     52 var ONE_DIV_SQRT_PI = 1.0 / SQRT_PI;
     53 var TWO_DIV_PI = 2.0 / PI;
     54 
     55 var x1 = 2.1971413260310170351e+00;
     56 var x2 = 5.4296810407941351328e+00;
     57 var x11 = 5.620e+02;
     58 var x12 = 1.8288260310170351490e-03;
     59 var x21 = 1.3900e+03;
     60 var x22 = -6.4592058648672279948e-06;
     61 
     62 // `sincos` workspace:
     63 var sc = [ 0.0, 0.0 ]; // WARNING: not thread safe
     64 
     65 
     66 // MAIN //
     67 
     68 /**
     69 * Computes the Bessel function of the second kind of order one.
     70 *
     71 * ## Notes
     72 *
     73 * -   Accuracy for subnormal `x` is very poor. Full accuracy is achieved at `1.0e-308` but trends progressively to zero at `5e-324`. This suggests that underflow (or overflow, perhaps due to a reciprocal) is effectively cutting off digits of precision until the computation loses all accuracy at `5e-324`.
     74 *
     75 * @param {number} x - input value
     76 * @returns {number} evaluated Bessel function
     77 *
     78 * @example
     79 * var v = y1( 0.0 );
     80 * // returns -Infinity
     81 *
     82 * v = y1( 1.0 );
     83 * // returns ~-0.781
     84 *
     85 * v = y1( -1.0 );
     86 * // returns NaN
     87 *
     88 * v = y1( Infinity );
     89 * // returns 0.0
     90 *
     91 * v = y1( -Infinity );
     92 * // returns NaN
     93 *
     94 * v = y1( NaN );
     95 * // returns NaN
     96 */
     97 function y1( x ) {
     98 	var rc;
     99 	var rs;
    100 	var y2;
    101 	var r;
    102 	var y;
    103 	var z;
    104 	var f;
    105 
    106 	if ( x < 0.0 ) {
    107 		return NaN;
    108 	}
    109 	if ( x === 0.0 ) {
    110 		return NINF;
    111 	}
    112 	if ( x === PINF ) {
    113 		return 0.0;
    114 	}
    115 	if ( x <= 4.0 ) {
    116 		y = x * x;
    117 		z = ( ln( x/x1 ) * besselj1( x ) ) * TWO_DIV_PI;
    118 		r = poly1( y );
    119 		f = ( ( x+x1 ) * ( (x - (x11/256.0)) - x12 ) ) / x;
    120 		return z + ( f*r );
    121 	}
    122 	if ( x <= 8.0 ) {
    123 		y = x * x;
    124 		z = ( ln( x/x2 ) * besselj1( x ) ) * TWO_DIV_PI;
    125 		r = poly2( y );
    126 		f = ( ( x+x2 ) * ( (x - (x21/256.0)) - x22 ) ) / x;
    127 		return z + ( f*r );
    128 	}
    129 	y = 8.0 / x;
    130 	y2 = y * y;
    131 	rc = polyC( y2 );
    132 	rs = polyS( y2 );
    133 	f = ONE_DIV_SQRT_PI / sqrt( x );
    134 
    135 	/*
    136 	* This code is really just:
    137 	*
    138 	* ```
    139 	* z = x - 0.75 * PI;
    140 	* return f * (rc * sin(z) + y * rs * cos(z));
    141 	* ```
    142 	*
    143 	* But using the sin/cos addition rules, plus constants for sin/cos of `3π/4` which then cancel out with corresponding terms in "f".
    144 	*/
    145 	sincos( sc, x );
    146 	return f * ( ( ( (y*rs) * (sc[0]-sc[1]) ) - ( rc * (sc[0]+sc[1]) ) ) );
    147 }
    148 
    149 
    150 // EXPORTS //
    151 
    152 module.exports = y1;