It would be good to get @andreasnoackjensen's feedback on this.

Awesome, thanks!

@johnmyleswhite I can confirm the slight difference between as121 and R/MATLAB/Mathematica. R cites: Amos, D. E. (1983). A portable Fortran subroutine for derivatives of the psi function, Algorithm 610, and it is part of SLATEC, http://www.netlib.org/slatec/src/dpsifn.f and therefore in the public domain. I tried to translate that to Julia at it seem to work. I think it is the more precise of the two algorithms and at the same time it gives all the derivatives so I suggest we use Amos' algorithm.

I don't have much to say about the inverse digamma function. I haven't encountered the function before and it also seems a little rare with around 240 hits on Google, but the implementation seems to work fine.

Do you have a working implementation of Amos' algorithm? If so, let's definitely use that.

The inverse digamma function comes up in fixed point algorithms involving the digamma function, e.g. MLE for the Dirichlet distribution.

I think it works. I'll make a pull request in a minute. There are still some open questions regarding the code, but it is easier to discuss when you can see the code.

These algorithms (especially invdigamma) look insanely slow compared to just chopping up the domain into a few pieces and using a minimax rational approximant on each piece. (Nowadays, Maple or Mathematica can construct minimax approximants for you, to arbitrary precision, so implementing this strategy is not so hard.)

But I supposed they can always be replaced by an optimized implementation later.

Closing this. I'll add invdigamma on its own now that we have polygamma.