Add Julia translation of psifn from SLATEC and some unicode constants

[andreasnoack]

This request is related to #3226

Right now I use i1mach and d1mach from opemlibm-extras since they were used in the SLATEC original but is there a more Julian way of getting those numbers? Also, any thoughts on the error messages. Right now most of them are taken from the source.

I allowed myself to introduce the Euler gamma for normal floats because it was already exported for BogFloats. I also introduced the unicode aliases π and γ. It has been discussed at some point and I like it.

[JeffBezanson]

It would be nice (but not essential) to translate the d1mach and i1mach calls into something more meaningful, For example eps, realmin, realmax, etc.

[johnmyleswhite]

This all looks good to me. I'd say merge this. After that I can make a new request for invdigamma. After that we can try to implement the strategy suggested by SGJ.

[johnmyleswhite]

One quick request: can we call this psigamma like R does? psifn makes me wish it was called psi alone, which is too useful a name to be stealing.

[StefanKarpinski]

-1 on "psigamma". Googling "trigamma" returns the wikipedia page about the appropriate function as the top hit. Googling "psigamma" returns a bunch of links to pages about sororities. Just because R chose a lousy, non-standard name doesn't mean we should make the same mistake.

[johnmyleswhite]

How about polygamma (http://www.mathworks.com/help/matlab/ref/psi.html)?

[StefanKarpinski]

Polygamma is the generalized function. The first, second, third and fourth specific instances are called "gamma", "digamma", "trigamma" and "tetragamma". It would seem weird to call the third one "psigamma", no?

[johnmyleswhite]

Andreas has implemented the generalized function and defined digamma and trigamma in terms of it. See lines 750 and 755.

[StefanKarpinski]

That seems quite sane to me. I guess we could always just have polygamma and say that you call that with appropriate n. The main advantage of having specific names would be if there are better/faster implementations of the special versions which can be called directly and which the generalized version can call for those n.

[johnmyleswhite]

I think the argument for having digamma and trigamma as separate is that most statistical applications will only need the first and second derivatives of the gamma function for defining things like the gradient of a distribution or its Hessian. You'll see that there are many real Google results for digamma and trigamma, but that tetragamma is quite sparse and there's really nobody talking about pentagamma.

[StefanKarpinski]

Agreed – we should definitely not have tetragamma, I was just pointing out that calling one of that sequence of functions psigamma is just strange.

[andreasnoack]

@JeffBezanson The function uses i1mach(14), i1mach(15) and i1mach(16) which are "t, the number of base-b digits" and "the smallest and largest exponent e". It also uses d1mach(4) and d1mach(5) which are "b**(1-t), the largest relative spacing" and "log10(b)". What would that be in Julia?

I like polygamma but there is a slight problem because psifn returns ((-1)^(k+1)/Γ(k+1))*ψ(k,x) where ψ(k,x) is the polygamma function. This is also a reason for having the di- and trigamma functions separately.

[StefanKarpinski]

I like polygamma but there is a slight problem because psifn returns ((-1)^(k+1)/Γ(k+1))*ψ(k,x) where ψ(k,x) is the polygamma function. This is also a reason for having the di- and trigamma functions separately.

I'm kind of confused by this statement. Can you elaborate?

[andreasnoack]

If Julia had a function polygamma(k,x) I guess people would expect it to return the (k+1)'th derivative of log(Γ(x)). The function psifn returns the scaled version mentioned earlier. We could of course just set polygamma(k,x)=(-1)^(k+1)*psifn(x,k,1,1)[1]/gamma(k+1) but I guess there is a good reason for defining psifn the way it is.

[johnmyleswhite]

I think exposing polygamma(k, x) as a method for accessing the unnormalized psifn would be the best strategy since people aren't going to know what to do with the last two arguments anyway. (I barely do.)

[andreasnoack]

I have hidden psifn and made an interface function polygamma which behaves as one would expect.

[johnmyleswhite]

I say go ahead with this. Trying to write a faster implementation later will be a fun exercise.

[johnmyleswhite]

I'm going to merge this now since it's holding up work on Distributions. The Travis build failure isn't related to math.

[johnmyleswhite]

@andreasnoackjensen, do you want to delete this branch now that it's merged?

[johnmyleswhite]

Just saw you did. Sorry for the spam.