<div dir="ltr"><div>Well crap ... I think I searched the Racket doc for digamma found nada, tried psi and blew it as I never progressed beyond "implements a polygamma function" and failed to recognize it for what it was. In other words, I flat missed it. So I was sort of asking about the process of voir dire'ing my impl with the idea of addition to the standard Racket math libraries, thinking it was absent, but ... never mind :). </div><div>Worthwhile I think to get 'digamma' to show up as a hit in the Racket doc search and maybe extending the note in the psi section, "Note that (psi 0 x) = (psi0 x)" to maybe "Note that (<a class="">psi</a> 0 x) <a class="">=</a> (<a class="">psi0</a> x) is the digamma function" for when someone searches on 'psi'. </div><div> </div></div><div class="gmail_extra"> <div class="gmail_quote">On Tue, May 7, 2013 at 11:39 AM, Neil Toronto <<a href="mailto:neil.toronto@gmail.com" target="_blank">neil.toronto@gmail.com</a>> wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im">On 05/06/2013 01:25 PM, Ray Racine wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> Neil, Here in is an implementation of the digamma function. </blockquote> </div> Cool! Welcome to the High Priesthood of the... crazy people who do floating-point implementations. :) BTW, I'm not clear on what you're requesting from me. Nearest I can come up with is how to measure error, so I'll address that below.<div class="im"> <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> Initially I tried Haskell's, then a Java impl from Mallet (LDA library) but they didn't align (BTW I think the Haskell impl is flat wrong for digamma (2)) and sought alternate sources. I finally found <a href="http://code.google.com/p/pmtk3/source/browse/trunk/util/digamma.m?r=227&spec=svn227" target="_blank">http://code.google.com/p/pmtk3/source/browse/trunk/util/digamma.m?r=227&spec=svn227</a> from pmtk2 a probabilistic modeling toolkit for Matlab/Octave, which is MIT LICENSED so I think we are good there. </blockquote> </div> You've discovered two things I found very frustrating while looking for code to copy for `math/special-functions' and `math/distributions': 1. Most implementations have ambiguous or terrible license terms. 2. Almost all implementations aren't good. #1 made me write a lot of new code. #2 let me feel okay with that. :D<div class="im"> <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> My impl is based (well ok copied) from that impl. Octave also appears to have a Gnu Science Library FFI which has an impl for psi (digamma) as well which I've used as a sanity check for positive x values. </blockquote> </div> You can use `bfpsi0' from `math/bigfloat' to get arbitrarily accurate error measurements. I used it to produce error plots when implementing `flpsi0'. Demo below.<div class="im"> <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> Since it works for x = 10 to a few decimal places, I'm sure it works everywhere ... perfectly. </blockquote> </div> :D<div class="im"> <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> So far the Racket impl appears to be aligned to Octave+GSL, but I have no intuition as what is acceptable epsilon variance in floating point vs what is true error in implementation. The impl handles z < 0, yet Octave+GSL blows up. </blockquote> </div> I set some pretty high standards for the math library. I call a flonum function implementation acceptable when: 1. It gives sensible answers on the largest possible domain. (So Octave+GSL's digamma doesn't qualify.) 2. Its error in ulps is <= 10 over most of the domain. 3. When error in ulps is higher in some subdomain, it's because there are no reasonably fast algorithms for computing the function accurately in that subdomain. Look up `flulp' and `flulp-error' in the math docs for definitions. Shortly, an ulp is the magnitude of the smallest change you can make to a flonum. Error in ulps is more or less the number of flonums between the true answer and the approximation. A function with error <= 0.5 ulps is *correct*. I call a function with <= 2 ulps friggin' fantastic, and one with <= 10 accurate. <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im"> I added addition poly coefficients up to 12 but the impl currently only uses 7. </div> impl follows: [...] </blockquote> Paste this at the bottom to test it: (require plot/typed math) (: digamma* (Flonum -> Real)) (define (digamma* x) (parameterize ([bf-precision 53]) (bigfloat->real (bfpsi0 (bf x))))) (: plot-error ((Flonum -> Flonum) (Flonum -> Real) Real Real -> Any)) (define (plot-error f-appx f-truth x-min x-max) (plot (function (λ: ([x : Real]) (let ([x (fl x)]) (define err (flulp-error (f-appx x) (f-truth x))) (if (rational? err) err -1.0))) x-min x-max))) (plot-error flpsi0 digamma* 0 32) (plot-error digamma digamma* 0 32) (plot-error flpsi0 digamma* -8 0) (plot-error digamma digamma* -8 0) It looks like your `digamma' is as accurate as `flpsi0' for x > 21 or so (i.e. <= 2 ulps, which is friggin' fantastic). If using 12 terms instead of 7 in your polynomial approximation reduced the error for smaller x, it would be worth using that polynomial approximation in the math library. (I can't test this because I don't know the coefficients' signs.) In the subdomain it would likely be accurate in, `flpsi0' uses 12- to 15-degree Chebyshev polynomials, which are about twice as slow to evaluate as similar-degree polynomials using Horner's rule (as in your `digamma-poly'). Neil ⊥ ____________________ Racket Users list: <a href="http://lists.racket-lang.org/users" target="_blank">http://lists.racket-lang.org/users</a> </blockquote></div> </div>