digamma.js (5134B)
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 http://www.boost.org/doc/libs/1_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/digamma.html}. The implementation follows the original but has been modified for JavaScript. 22 * 23 * ```text 24 * (C) Copyright John Maddock 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 isnan = require( './../../../../base/assert/is-nan' ); 37 var floor = require( './../../../../base/special/floor' ); 38 var tan = require( './../../../../base/special/tan' ); 39 var PI = require( '@stdlib/constants/float64/pi' ); 40 var asymptoticApprox = require( './asymptotic_expansion.js' ); 41 var rationalApprox = require( './rational_approximation.js' ); 42 43 44 // VARIABLES // 45 46 var MIN_SAFE_ASYMPTOTIC = 10.0; // BIG! 47 48 49 // MAIN // 50 51 /** 52 * Evaluates the digamma function. 53 * 54 * ## Method 55 * 56 * 1. For \\(x < 0\\), we use the reflection formula 57 * 58 * ```tex 59 * \psi(1-x) = \psi(x) + \frac{\pi}{\tan(\pi x)} 60 * ``` 61 * 62 * to make \\(x\\) positive. 63 * 64 * 2. For \\(x \in \[0,1]\\), we use the recurrence relation 65 * 66 * ```tex 67 * \psi(x) = \psi(x+1) - \frac{1}{x} 68 * ``` 69 * 70 * to shift the evaluation range to \\(\[1,2]\\). 71 * 72 * 3. For \\(x \in \[1,2]\\), we use a rational approximation of the form 73 * 74 * ```tex 75 * \psi(x) = (x - \mathrm{root})(Y + \operatorname{R}(x-1)) 76 * ``` 77 * 78 * where \\(\mathrm{root}\\) is the location of the positive root of \\(\psi\\), \\(Y\\) is a constant, and \\(R\\) is optimized for low absolute error compared to \\(Y\\). 79 * 80 * <!-- <note>--> 81 * 82 * Note that, since \\(\mathrm{root}\\) is irrational, we need twice as many digits in \\(\mathrm{root}\\) as in \\(x\\) in order to avoid cancellation error during subtraction, assuming \\(x\\) has an exact value. This means that, even if \\(x\\) is rounded to the next representable value, the result of \\(\psi(x)\\) will not be zero. 83 * 84 * <!-- </note> --> 85 * 86 * <!-- <note> --> 87 * 88 * This approach gives 17-digit precision. 89 * 90 * <!-- </note> --> 91 * 92 * 4. For \\(x \in \[2,\mathrm{BIG}]\\), we use the recurrence relation 93 * 94 * ```tex 95 * \psi(x+1) = \psi(x) + \frac{1}{x} 96 * ``` 97 * 98 * to shift the evaluation range to \\(\[1,2]\\). 99 * 100 * 5. For \\(x > \mathrm{BIG}\\), we use the asymptotic expression 101 * 102 * ```tex 103 * \psi(x) = \ln(x) + \frac{1}{2x} - \biggl( \frac{B_{21}}{2x^2} + \frac{B_{22}}{4x^4} + \frac{B_{23}}{6x^6} + \ldots \biggr) 104 * ``` 105 * 106 * This expansion, however, is divergent after a few terms. The number of terms depends on \\(x\\). Accordingly, we must choose a value of \\(\mathrm{BIG}\\) which allows us to truncate the series at a term that is too small to have an effect on the result. Setting \\(\mathrm{BIG} = 10\\), allows us to truncate the series early and evaluate as \\(1/x^2\\). 107 * 108 * <!-- <note> --> 109 * 110 * This approach gives 17-digit precision for \\(x \geq 10\\). 111 * 112 * <!-- </note> --> 113 * 114 * ## Notes 115 * 116 * - Maximum deviation found: \\(1.466\\mbox{e-}18\\) 117 * - Max error found: \\(2.452\mbox{e-}17\\) (double precision) 118 * 119 * 120 * @param {number} x - input value 121 * @returns {number} function value 122 * 123 * @example 124 * var v = digamma( -2.5 ); 125 * // returns ~1.103 126 * 127 * @example 128 * var v = digamma( 1.0 ); 129 * // returns ~-0.577 130 * 131 * @example 132 * var v = digamma( 10.0 ); 133 * // returns ~2.252 134 * 135 * @example 136 * var v = digamma( NaN ); 137 * // returns NaN 138 * 139 * @example 140 * var v = digamma( -1.0 ); 141 * // returns NaN 142 */ 143 function digamma( x ) { 144 var rem; 145 var tmp; 146 if ( isnan( x ) || x === 0.0 ) { 147 return NaN; 148 } 149 // If `x` is negative, use reflection... 150 if ( x <= -1.0 ) { 151 // Reflect: 152 x = 1.0 - x; 153 154 // Argument reduction for tan: 155 rem = x - floor(x); 156 157 // Shift to negative if > 0.5: 158 if ( rem > 0.5 ) { 159 rem -= 1.0; 160 } 161 // Check for evaluation at a negative pole: 162 if ( rem === 0.0 ) { 163 return NaN; 164 } 165 tmp = PI / tan( PI * rem ); 166 } else { 167 tmp = 0.0; 168 } 169 // If we're above the lower-limit for the asymptotic expansion, then use it... 170 if ( x >= MIN_SAFE_ASYMPTOTIC ) { 171 tmp += asymptoticApprox( x ); 172 return tmp; 173 } 174 // If x > 2, reduce to the interval [1,2]... 175 while ( x > 2.0 ) { 176 x -= 1.0; 177 tmp += 1.0/x; 178 } 179 // If x < 1, use recurrence to shift to > 1.. 180 while ( x < 1.0 ) { 181 tmp -= 1.0/x; 182 x += 1.0; 183 } 184 tmp += rationalApprox( x ); 185 return tmp; 186 } 187 188 189 // EXPORTS // 190 191 module.exports = digamma;