Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 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 125 126 127 128 129 130 131 132 133 134 | 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 917x 805x 383x 383x 805x 805x 112x 917x 320x 320x 320x 320x 320x 320x 112x 112x 208x 208x 208x 208x 208x 320x 320x 112x 917x 112x 917x 917x 1x 1x 1x 1x 1x | /**
* @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 original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_65_0/boost/math/special_functions/detail/polygamma.hpp}. The implementation follows the original but has been modified for JavaScript.
*
* ```text
* (C) Copyright Nikhar Agrawal 2013.
* (C) Copyright Christopher Kormanyos 2013.
* (C) Copyright John Maddock 2014.
* (C) Copyright Paul Bristow 2013.
*
* Use, modification and distribution are subject to the
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
* ```
*/
'use strict';
// MODULES //
var logger = require( 'debug' );
var factorial = require( '@stdlib/math/base/special/factorial' );
var zeta = require( '@stdlib/math/base/special/riemann-zeta' );
var abs = require( '@stdlib/math/base/special/abs' );
var pow = require( '@stdlib/math/base/special/pow' );
var NINF = require( '@stdlib/constants/float64/ninf' );
var PINF = require( '@stdlib/constants/float64/pinf' );
var EPS = require( '@stdlib/constants/float64/eps' );
var MAX = require( '@stdlib/constants/float64/max' );
// VARIABLES //
var debug = logger( 'polygamma' );
var MAX_SERIES_ITERATIONS = 1000000;
// MAIN //
/**
* Evaluates the polygamma function near zero.
*
* ## Notes
*
* - If we take this [expansion][1] for `polygamma` and substitute in this [expression][2] for `polygamma(n, 1)`, we get an alternating series for polygamma when `x` is small in terms of zeta functions of integer arguments (which are easy to evaluate, at least when the integer is even).
*
* [1]: http://functions.wolfram.com/06.15.06.0003.02
* [2]: http://functions.wolfram.com/06.15.03.0009.01
*
* @private
* @param {PositiveInteger} n - derivative to evaluate
* @param {number} x - input value
* @returns {number} (n+1)'th derivative
*/
function nearzero( n, x ) {
var factorialPart;
var prefix;
var scale;
var term;
var sum;
var AX;
var k;
// In order to avoid spurious overflow, save the `n!` term for later, and rescale at the end:
scale = factorial( n );
// "factorialPart" contains everything except the zeta function evaluations in each term:
factorialPart = 1;
// "prefix" is what we'll be adding the accumulated sum to, it will be `n! / z^(n+1)`, but since we're scaling by `n!` it is just `1 / z^(n+1)` for now:
prefix = pow( x, n+1 );
if ( prefix === 0.0 ) {
return PINF;
}
prefix = 1.0 / prefix;
// First term in the series is necessarily `< zeta(2) < 2`, so ignore the sum if it will have no effect on the result:
if ( prefix > 2.0/EPS ) {
if ( n & 1 ) {
return ( AX/prefix < scale ) ? PINF : prefix * scale;
}
return ( AX/prefix < scale ) ? NINF : -prefix * scale;
}
sum = prefix;
for ( k = 0; ; ) {
// Get the k'th term:
term = factorialPart * zeta( k+n+1 );
sum += term;
// Termination condition:
if ( abs( term ) < abs(sum * EPS ) ) {
break;
}
// Move on `k` and `factorialPart`:
k += 1;
factorialPart *= (-x * (n+k)) / k;
// Last chance exit:
if ( k > MAX_SERIES_ITERATIONS ) {
debug( 'Series did not converge, best value is %d.', sum );
return NaN;
}
}
// We need to multiply by the scale, at each stage checking for overflow:
if ( MAX/scale < sum ) {
return PINF;
}
sum *= scale;
return ( n & 1 ) ? sum : -sum;
}
// EXPORTS //
module.exports = nearzero;
|