public static double s14ac(
	double x,
	out int ifail
)

Public Shared Function s14ac ( _
	x As Double, _
	<OutAttribute> ByRef ifail As Integer _
) As Double

public:
static double s14ac(
	double x, 
	[OutAttribute] int% ifail
)

static member s14ac : 
        x : float * 
        ifail : int byref -> float 

s14ac returns a value of the function $\psi \left(x\right)-\ln x$, where $\psi$ is the psi function $\psi \left(x\right)=\frac{d}{dx}\ln \Gamma \left(x\right)=\frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}$.

s14ac returns a value of the function $\psi \left(x\right)-\ln x$. The psi function is computed without the logarithmic term so that when $x$ is large, sums or differences of psi functions may be computed without unnecessary loss of precision, by analytically combining the logarithmic terms. For example, the difference $d=\psi \left(x+\frac{1}{2}\right)-\psi \left(x\right)$ has an asymptotic behaviour for large $x$ given by $d\sim \ln \left(x+\frac{1}{2}\right)-\ln x+\mathit{O}\left(\frac{1}{x^2}\right)\sim \ln \left(1+\frac{1}{2x}\right)\sim \frac{1}{2x}$.

Computing $d$ directly would amount to subtracting two large numbers which are close to $\ln \left(x+\frac{1}{2}\right)$ and $\ln x$ to produce a small number close to $\frac{1}{2x}$, resulting in a loss of significant digits. However, using this method to compute $f\left(x\right)=\psi \left(x\right)-\ln x$, we can compute $d=f\left(x+\frac{1}{2}\right)-f\left(x\right)+\ln \left(1+\frac{1}{2x}\right)$, and the dominant logarithmic term may be computed accurately from its power series when $x$ is large. Thus we avoid the unnecessary loss of precision.

The method is derived from the method PSIFN in Amos (1983).

Abramowitz M and Stegun I A (1972) Handbook of Mathematical Functions (3rd Edition) Dover Publications

Amos D E (1983) Algorithm 610: A portable FORTRAN subroutine for derivatives of the psi function ACM Trans. Math. Software 9 494–502

Errors or warnings detected by the method:

All constants in s14ac are given to approximately $18$ digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used $t$, then clearly the maximum number of correct digits in the results obtained is limited by $p=\min \left(t,18\right)$.

With the above proviso, results returned by this method should be accurate almost to full precision, except at points close to the zero of $\psi \left(x\right)$, $x\simeq 1.461632$, where only absolute rather than relative accuracy can be obtained.

None.

None.

The example program reads values of the argument $x$ from a file, evaluates the function at each value of $x$ and prints the results.