R devel - Feb 2004 - digamma with negative arguments (PR#6626)

Full_Name: Chanseok Park
Version: 1.8.1
OS: linux-gnu
Submission from: (NULL) (130.127.112.183)


digamma with any negative value does not give a right answer.
It gives -1.797693e+308 for any negative arguments.

For example, digamma(-1.1) gives  -1.797693e+308.

The right answer should be 10.15416

This bug can be easily fixed by using the following digamma identity.

digamma(x) = digamma(1-x) - pi/tan(pi*x) ## if x is negative
> digamma(-1)
[1] -1.797693e+308

## if we use the above math identity, we have the right answer 
##     for any negative arguments.
> 
> x= -1.1
>  digamma(1-x) - pi/tan(pi*x) ## if x is negative
[1] 10.15416
>

I intend to look into this :
>>>>> "cspark" == cspark 
<cspark@ces.clemson.edu>
>>>>>     on Sun, 29 Feb 2004 22:15:02 +0100 (CET) writes:
    cspark> Full_Name: Chanseok Park

    cspark> digamma with any negative value does not give a right answer.
    cspark> It gives -1.797693e+308 for any negative arguments.

    cspark> For example, digamma(-1.1) gives  -1.797693e+308.

    cspark> The right answer should be 10.15416

    cspark> This bug can be easily fixed by using the following digamma
identity.
    cspark> digamma(x) = digamma(1-x) - pi/tan(pi*x) ## if x is negative

    >> digamma(-1)
    cspark> [1] -1.797693e+308

    cspark> ## if we use the above math identity, we have the right answer 
    cspark> ##     for any negative arguments.
    >> 
    >> x= -1.1
    >> digamma(1-x) - pi/tan(pi*x) ## if x is negative
    cspark> [1] 10.15416
    >> 

and there are similar "reflection formulae" for trigamma() etc.
which I'd want as well.

As a matter of fact, our internal code src/nmath/polygamma.c
uses the basic workhorse  dpsifn()  which provides "arbitrary"
derivatives of the psi() { == digamma() } function {for positive x}.

dpsifn() needs input/output arrays in general,

  static void
  dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr)

but I think it would still be worth to remove "static", add
documentation in R-ext.texi and hence add it to the API.
Another step would be an R-level interface; maybe with
   digamma(x, deriv = 0)
	      ^^^^^^^^^
with obvious meaning and implementation.

Feedback ?
{divert from R-bugs to R-devel, please !}

Martin Maechler <maechler@stat.math.ethz.ch>
http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><