Re: psi(1,x) (was: Making some minor contributions to GSL)

[Gerard Jungman <jungman at lanl dot gov>]

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