<div dir="ltr"><div>Well crap ... <br><br>I think I searched the Racket doc for digamma found nada, tried psi and blew it as I never progressed beyond &quot;implements a polygamma function&quot; 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&#39;ing my impl with the idea of addition to the standard Racket math libraries, thinking it was absent, but ... never mind :).   <br>
<br></div><div>Worthwhile I think to get &#39;digamma&#39; to show up as a hit in the Racket doc search and maybe extending the note in the psi section, &quot;Note that (psi 0 x) = (psi0 x)&quot; to maybe &quot;Note that
<span class="">(</span><span class=""><a class="">psi</a></span><span class=""> </span><span class="">0</span><span class=""> </span><span class="">x</span><span class="">)</span><span class=""> </span><span class=""><a class="">=</a></span><span class=""> </span><span class="">(</span><span class=""><a class="">psi0</a></span><span class=""> </span><span class="">x</span><span class="">)</span> is the digamma function&quot; for when someone searches on &#39;psi&#39;.<br>
<br></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, May 7, 2013 at 11:39 AM, Neil Toronto <span dir="ltr">&lt;<a href="mailto:neil.toronto@gmail.com" target="_blank">neil.toronto@gmail.com</a>&gt;</span> wrote:<br>
<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:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Neil,<br>
<br>
Here in is an implementation of the digamma function.<br>
</blockquote>
<br></div>
Cool! Welcome to the High Priesthood of the... crazy people who do floating-point implementations. :)<br>
<br>
BTW, I&#39;m not clear on what you&#39;re requesting from me. Nearest I can come up with is how to measure error, so I&#39;ll address that below.<div class="im"><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Initially I tried<br>
Haskell&#39;s, then a Java impl from Mallet  (LDA library) but they didn&#39;t<br>
align (BTW I think the Haskell impl is flat wrong for digamma (2)) and<br>
sought alternate sources.<br>
<br>
I finally found<br>
<a href="http://code.google.com/p/pmtk3/source/browse/trunk/util/digamma.m?r=227&amp;spec=svn227" target="_blank">http://code.google.com/p/<u></u>pmtk3/source/browse/trunk/<u></u>util/digamma.m?r=227&amp;spec=<u></u>svn227</a><br>

from pmtk2 a probabilistic modeling toolkit for Matlab/Octave, which is<br>
MIT LICENSED so I think we are good there.<br>
</blockquote>
<br></div>
You&#39;ve discovered two things I found very frustrating while looking for code to copy for `math/special-functions&#39; and `math/distributions&#39;:<br>
<br>
1. Most implementations have ambiguous or terrible license terms.<br>
<br>
2. Almost all implementations aren&#39;t good.<br>
<br>
#1 made me write a lot of new code. #2 let me feel okay with that. :D<div class="im"><br>
<br>
<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.<br>
<br>
Octave also appears to have a Gnu Science Library FFI which has an impl<br>
for psi (digamma) as well which I&#39;ve used as a sanity check for positive<br>
x values.<br>
</blockquote>
<br></div>
You can use `bfpsi0&#39; from `math/bigfloat&#39; to get arbitrarily accurate error measurements. I used it to produce error plots when implementing `flpsi0&#39;. Demo below.<div class="im"><br>
<br>
<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&#39;m sure it works<br>
everywhere ... perfectly.<br>
</blockquote>
<br></div>
:D<div class="im"><br>
<br>
<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<br>
no intuition as what is acceptable epsilon variance in floating point vs<br>
what is true error in implementation.  The impl handles z &lt; 0, yet<br>
Octave+GSL blows up.<br>
</blockquote>
<br></div>
I set some pretty high standards for the math library. I call a flonum function implementation acceptable when:<br>
<br>
 1. It gives sensible answers on the largest possible domain. (So<br>
    Octave+GSL&#39;s digamma doesn&#39;t qualify.)<br>
<br>
 2. Its error in ulps is &lt;= 10 over most of the domain.<br>
<br>
 3. When error in ulps is higher in some subdomain, it&#39;s because there<br>
    are no reasonably fast algorithms for computing the function<br>
    accurately in that subdomain.<br>
<br>
Look up `flulp&#39; and `flulp-error&#39; 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.<br>

<br>
A function with error &lt;= 0.5 ulps is *correct*. I call a function with &lt;= 2 ulps friggin&#39; fantastic, and one with &lt;= 10 accurate.<br>
<br>
<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<br>
uses 7.<br>
<br></div>
impl follows:  [...]<br>
</blockquote>
<br>
Paste this at the bottom to test it:<br>
<br>
<br>
(require plot/typed<br>
         math)<br>
<br>
(: digamma* (Flonum -&gt; Real))<br>
(define (digamma* x)<br>
  (parameterize ([bf-precision  53])<br>
    (bigfloat-&gt;real (bfpsi0 (bf x)))))<br>
<br>
(: plot-error ((Flonum -&gt; Flonum) (Flonum -&gt; Real) Real Real -&gt; Any))<br>
(define (plot-error f-appx f-truth x-min x-max)<br>
  (plot (function<br>
         (λ: ([x : Real])<br>
           (let ([x  (fl x)])<br>
             (define err (flulp-error (f-appx x) (f-truth x)))<br>
             (if (rational? err) err -1.0)))<br>
         x-min x-max)))<br>
<br>
(plot-error flpsi0 digamma* 0 32)<br>
(plot-error digamma digamma* 0 32)<br>
<br>
(plot-error flpsi0 digamma* -8 0)<br>
(plot-error digamma digamma* -8 0)<br>
<br>
<br>
It looks like your `digamma&#39; is as accurate as `flpsi0&#39; for x &gt; 21 or so (i.e. &lt;= 2 ulps, which is friggin&#39; 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&#39;t test this because I don&#39;t know the coefficients&#39; signs.) In the subdomain it would likely be accurate in, `flpsi0&#39; uses 12- to 15-degree Chebyshev polynomials, which are about twice as slow to evaluate as similar-degree polynomials using Horner&#39;s rule (as in your `digamma-poly&#39;).<br>

<br>
Neil ⊥<br>
<br>
____________________<br>
 Racket Users list:<br>
 <a href="http://lists.racket-lang.org/users" target="_blank">http://lists.racket-lang.org/<u></u>users</a><br>
</blockquote></div><br></div>