Numerical computation of Tricomi's psi function by the trapezoidal rule (Q1822445)

The trapezoidal rule is applied to the numerical calculation of the integral representation of Tricomi's psi function \(\Psi (a,c;x)=2x^{1- c}/\Gamma (a)\int ^{+\infty}_{0}e^{-u^ 2}(u^{2a-1}/(x+u^ 2)^{1-c+a})du\) for \(a,x\in {\mathbb{R}}^ +\), \(c\in {\mathbb{R}}\). The unexpectedly high accuracy is explained by means of a careful investigation in the complex field of the Euler-Maclaurin formula, and particularly of its remainder terms, considered as an extension of the trapezoidal rule. Since the same method has been used previously to evaluate the complementary incomplete gamma function [the authors, Numer. Math. 50, 419-428 (1987; Zbl 0593.65017)], the Euler gamma function and the digamma function, the present paper limits itself to quote the main features of the method and describes fully the peculiarities of this application. A simple and efficient numerical procedure for obtaining values of the psi function is given; moreover, to reduce the amount of calculation, an iterative algorithm for the evaluation of the trapezoidal rule, similar to Horner's scheme for polynomials, is suggested.