A little code study, formula study and experimentation reveals that the
situation is mostly fixable:
1. Get rid of the explicit alpha limit. (A form of it is implicit in
(2) and (3) below.)
2. Use the series formula when
(x < alph + 10 && x < 0.99 * (alph + 1000))
This guarantees that the sum converges reasonably fast. (For extremely
large x and alpha that will take about -53/log2(0.99) iterations for
53 significant bits, i.e., about 3700 iterations.)
3. Use the continued fraction formula when
(alph < x && alph - 100 < 0.99 * x)
Aka, you don't want to use the formula either near the critical point
where alpha/x ~ 1 unless the numbers are small.
4a. Go to a library and figure out how Temme does it for alpha near x,
both large. In this case the 0.99 from above could probably be
lowered a lot for faster convergence.
or
4b. Use the pnorm approximation. It seems to do a whole lot better for
alpha near x than it did for the 10:1 case I quoted.
Comments, please.
Morten
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
