function value = trigamma ( x )

%*****************************************************************************80
%
%% TRIGAMMA calculates trigamma(x) = d**2 log(Gamma(x)) / dx**2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    B Schneider,
%    Trigamma Function,
%    Algorithm AS 121,
%    Applied Statistics,
%    Volume 27, Number 1, page 97-99, 1978.
%
%  Parameters:
%
%    Input, X, the argument of the trigamma function.
%    0 < X.
%
%    Output, real VALUE, the value of the trigamma function at X.
%
  a = 0.0001;
  b = 5.0;
  b2 =  1.0 / 6.0;
  b4 = -1.0 / 30.0;
  b6 =  1.0 / 42.0;
  b8 = -1.0 / 30.0;
%
%  1): If X is not positive, fail.
%
  if ( x <= 0.0 )

    value = 0.0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIGAMMA - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'TRIGAMMA - Fatal error!' );
%
%  2): If X is smaller than A, use a small value approximation.
%
  elseif ( x <= a )

    value = 1.0 / x^2;
%
%  3): Otherwise, increase the argument to B <= ( X + I ).
%
  else

    z = x;
    value = 0.0;

    while ( z < b )
      value = value + 1.0 / z^2;
      z = z + 1.0;
    end
%
%  ...and then apply an asymptotic formula.
%
    y = 1.0 / z^2;

    value = value + 0.5 * ...
            y + ( 1.0 ...
          + y * ( b2 ...
          + y * ( b4 ...
          + y * ( b6 ...
          + y *   b8 )))) / z;

  end

  return
end