function value = digamma ( x ) %*****************************************************************************80 % %% DIGAMMA calculates the digamma or Psi function = d ( LOG ( GAMMA ( X ) ) ) / dX % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 07 September 2004 % % Author: % % MATLAB version by John Burkardt. % % Reference: % % J Bernardo, % Psi ( Digamma ) Function, % Algorithm AS 103, % Applied Statistics, % Volume 25, Number 3, pages 315-317, 1976. % % Parameters: % % Input, real X, the argument of the digamma function. % 0 < X. % % Output, real VALUE, the value of the digamma function at X. % % % Check the input. % if ( x <= 0.0 ) value = 0.0; return end % % Use approximation for small argument. % if ( x <= 0.000001 ) euler_mascheroni = 0.57721566490153286060; value = - euler_mascheroni - 1.0 / x + 1.6449340668482264365 * x; return end % % Reduce to DIGAMA(X + N). % value = 0.0; while ( x < 8.5 ) value = value - 1.0 / x; x = x + 1.0; end % % Use Stirling's (actually de Moivre's) expansion. % r = 1.0 / x; value = value + log ( x ) - 0.5 * r; r = r * r; value = value ... - r * ( 1.0 / 12.0 ... - r * ( 1.0 / 120.0 ... - r * ( 1.0 / 252.0 ... - r * ( 1.0 / 240.0 ... - r * ( 1.0 / 132.0 ) ) ) ) ); return end