<div dir="ltr">And thanks for the overview on FP implementation and evaluation, very instructive. Saving this email for reference down the road. </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>