psi

Psi (digamma) function.

SYNOPSIS:

double x, y, psi();
y = psi(x);

DESCRIPTION:

              d      -
   psi(x) = -- ln | (x)
              dx

is the logarithmic derivative of the gamma function.
For integer x:
                   n-1
                    -
psi(n) = -EUL +   > 1/k.
                    -
                   k=1

This formula is used for 0 < n <= 10. If x is negative, it is transformed to a positive argument by the reflection formula psi(1-x) = psi(x) + pi cot(pi x). For general positive x, the argument is made greater than 10 using the recurrence psi(x+1) = psi(x) + 1/x. Then this asymptotic expansion is applied:

                           inf.   B
                            -      2k
psi(x) = log(x) - 1/2x -   >    -------
                            -        2k
                           k=1   2k x

where the B2k are Bernoulli numbers.

ACCURACY:

Relative error (except absolute when |psi| < 1):

arithmetic   domain     # trials      peak         rms
    DEC       0,30         2500       1.7e-16     2.0e-17
    IEEE      0,30        30000       1.3e-15     1.4e-16
    IEEE      -30,0       40000       1.5e-15     2.2e-16

ERROR MESSAGES:

message         condition      value returned
singularity    x integer <=0      MAXNUM