The issue seems to be that the infinite sum is truncated too early when x
is in the extreme upper tail. An easily validated improvement to to
dnchisq.c would be to add an additional requirement in the upper tail
while condition, that the summation should continue while the additive
term remains too big a fraction of the saddlepoint density. Here is
replacement code that partially offsets the time required for addtional
terms in the sum by calculating x*ncp2 only once instead of repeating it
in every iteration of both loops.
x2 = x * ncp2
/* upper tail */
term = mid; df = dfmid; i = imax;
s = 0.5 - (dx+sqrt(dx*dx+4*ncp/x))*0.25 /* saddlepoint */
s2 = 1-2*s
termlimit = exp(ncp*s/s2 - log(s2)*df*0.5 - s*x -
log((df+2*ncp+4*ncp*s/s2)/(s2*s2)*pi) - log(2)*2**54) /* saddlepoint
density * 2^-53 */
do {
i++;
q = x2 / i / df;
df += 2;
term *= q;
sum += term;
} while (q >= 1 || term * q > (1-q)*eps || term > termlimit);
/* lower tail */
term = mid; df = dfmid; i = imax;
while ( i ){
df -= 2;
q = i * df / x2;
i--;
term *= q;
sum += term;
if (q < 1 && term * q <= (1-q)*eps) break;
}
Jerry
[[alternative HTML version deleted]]