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
|
/**
* @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( './../../../../base/special/factorial' );
var zeta = require( './../../../../base/special/riemann-zeta' );
var abs = require( './../../../../base/special/abs' );
var pow = require( './../../../../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;
|
