Making some minor contributions to GSL

linas@austin.ibm.com linas@austin.ibm.com
Wed Jan 14 00:57:00 GMT 2004


On Tue, Jan 13, 2004 at 03:56:00PM -0700, Gerard Jungman wrote:
> On Tue, 2004-01-13 at 11:17, linas@austin.ibm.com wrote:
> > 
> > I have some code that might be appropriate for inclusion with GSL,
> > and was curious to gauge interest, and understand the hurdles 
> > involved.  
> > 
> > Polygamma functions, Abramowitz & Stegun chapter 6.4 -- these
> > are derivatives of the gamma (factorial) function.   The algo
> > I have is fairly straightforward ... and is accurate to about
> > 1e-14 or so for "small" orders and arguments.  (The "multiplication
> > theorm" breaks down above order 30 or 40 or so, i.e. for arguments
> > greater than a 100 or 200 or so. I'm not sure why.)
> 
> I have never made a study of robustness of the GSL
> implementation, in gsl_sf_psi_n(). There are only

How embarassing.  When I looked for it, I didn't find it, 
so I wrote my own implementation ... 

> the test cases in specfunc/test_sf.c. Is there
> a specific deficiency in the GSL implementation?

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.  

You may find it reassuring that your gsl_sf_psi_n (1,x) agrees
with mine to better than 1e-14 across the board, from very close 
to the pole at x=0 (1e-6)  to at least 1e+6.  I beleive that
most of this error is from my own algo, which wasn't tuned.

It might make you feel better to also know that your algo is 
faster than what I came up with, a lot faster in some cases,
especially for large x.

> If you have some test cases that break it, please
> send them so I can check it out. If you have code
> that can plug a hole in the current implementation,
> then that would be greatly appreciated.

Well, now that I know it exists ... it seems to be quite fine.

=====
Another minor suggestion/comment: It would be nice to have 
a function that returns Riemann zeta -1, which provides better
accuracy for the intermediate/larger values values of the argument. 

Note also, (numerically) convergent sums involving zeta are more 
likely to use it in the form of (zeta - 1), which falls off very 
nicely.  At least, that's the way I like it ...

--linas



More information about the Gsl-discuss mailing list