psi(1,x) (was: Making some minor contributions to GSL)
Gerard Jungman
jungman@lanl.gov
Thu Jan 22 00:20:00 GMT 2004
On Fri, 2004-01-16 at 17:24, linas@austin.ibm.com wrote:
> On Thu, Jan 15, 2004 at 05:57:31PM -0700, Gerard Jungman wrote:
> > On Tue, 2004-01-13 at 17:57, linas@austin.ibm.com wrote:
> >
> > > The only "deficiency" is that the GSL implementation generates
> > > a "domain error" when gsl_sf_psi_n (1,x) is called with
> > > negative x. The psi's should be perfectly well defined
> > > for negative, non-integer x's.
> >
> > It might take me a while to get around to that. I don't think
> > I ever found a good way, or at least not one that I understood.
> I'm not sure what you mean. For small negative values,
> abram & steg 6.4.6 can be applied a few times and for larger
> values, 6.4.7. I'm not sure, I guess you must be worried about
> precision/rounding or performance or something like that?
> Or maybe large values of n, for which this isn't practical?
Now in CVS:
gsl_sf_psi_1_e(double x, gsl_sf_result * r);
gsl_sf_psi_1_e(double x);
works for all x != 0, -1, -2, ...
Yes, what I meant was that I didn't know a good way for general n and x.
But as you point out, for small fixed n it is easy. Since we already
had a psi(1,n) for positive integer argument, I decided to also
implement psi(1,x), for parallelism with the psi(0,.) implementations.
Also, gsl_sf_psi(int n, double x) will now work for negative x, for
the cases n=0,1. This feature is un-documented for now, since there
is no good way to explain that it won't work for n > 1.
Eventually I may be able extend this to all x for n > 1. Perhaps the
best way is to extend gsl_sf_hzeta(s,q) to handle q < 0, q != 0, -1,
-2,... This is how I compute the general psi(n,x) case currently,
for x > 0, using A+S 6.4.10. Since we only need integer s in this use of
hzeta(s,x), it may not be too hard. I should probably implement that
restriction for hzeta directly, i.e. hzeta_int(int s, double q).
Anyway, since you mentioned the case psi(1,.) explicitly in your
message, I thought I could at least do that for now.
Thanks.
--
Gerard Jungman <jungman@lanl.gov>
Los Alamos National Laboratory
More information about the Gsl-discuss
mailing list