No subject

From: ()
Date: Mon Dec 3 19:58:15 EST 2012

|#


;; Assumes for z < 0
;;  Use the reflection formula (Jeffrey 11.1.6):
;;  digamma(-x) =3D digamma(x+1) + pi*cot(pi*x)
;;   y =3D digamma(-x+1) + pi*cot(-pi*x);
;;  This is related to the identity
;;  digamma(-x) =3D digamma(x+1) - digamma(z) + digamma(1-z)
;;  where z is the fractional part of x
;;  For example:
;;  digamma(-3.1) =3D 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1)
;;                =3D digamma(4.1) - digamma(0.1) + digamma(1-0.1)
;;  Then we use
;;  digamma(1-z) - digamma(z) =3D pi*cot(pi*z)
;;  cot(x) =3D tan (pi/2 - x)

(provide:
 [digamma (Float -> Float)])

(require
 racket/match
 racket/pretty
 (only-in racket/flonum
          fltan
          fllog)
 (only-in racket/math
          pi))

(define: PI : Float pi)
(define: NEG-PI : Float (- PI))
(define: PI-OVER-TWO : Float (/ PI 2.0))
(define: PI-SQUARED-OVER-SIX : Float (/ (* PI PI) 6))
(define: EULER-MASCHERONI : Float -0.5772156649015328606065121)

(define: DIGAMMA-1 : Float EULER-MASCHERONI)     ;; digamma(1)
(define: DIGAMMA-2 : Float PI-SQUARED-OVER-SIX)

(define: DIGAMMA-COEF-2 : Float (- (/ 1.0 12.0)))
(define: DIGAMMA-COEF-3 : Float (/ 1.0 120.0))
(define: DIGAMMA-COEF-4 : Float (- (/ 1.0 252.0)))
(define: DIGAMMA-COEF-5 : Float (/ 1.0 240.0))
(define: DIGAMMA-COEF-6 : Float (- (/ 1.0 132.0)))
(define: DIGAMMA-COEF-7 : Float (/ 691.0 32760.0))
(define: DIGAMMA-COEF-8 : Float (- (/ 1.0 12.0)))

(define: DIGAMMA-LARGE : Float 9.5)
(define: DIGAMMA-SMALL : Float 0.000001)

;;anxn + =E2=80=A6 + a2x2 + a1x + a0 =3D a0 + x(a1 + x(a2 + x(=E2=80=A6 an)=
 =E2=80=A6) )
(: digamma-poly (Float -> Float))
(define (digamma-poly z)
  (* z (+ DIGAMMA-COEF-2
          (* z (+ DIGAMMA-COEF-3
                  (* z (+ DIGAMMA-COEF-4
                          (* z (+ DIGAMMA-COEF-5
                                  (* z (+ DIGAMMA-COEF-6
                                          (* z (+ DIGAMMA-COEF-7
                                                  (* z (+
DIGAMMA-COEF-8)))))))))))))))

(: cotan (Float -> Float))
(define (cotan x)
  (fltan (- PI-OVER-TWO x)))

(: digamma-negative (Float -> Float))
(define (digamma-negative z)
  ;; (assert (< z 0.0))
  (+ (digamma (+ (- z) 1.0))
     (* PI (cotan (* z NEG-PI)))))

(: digamma-shift (Float -> (Values Float Float)))
(define (digamma-shift z)
  (let loop ((z z) (y 0.0))
    (if (>=3D z DIGAMMA-LARGE)
        (values z y)
        (loop (+ z 1.0) (- y (/ z))))))

(: digamma (Float -> Float))
(define (digamma z)
  (cond
    ((eq? z 0.0)
     -inf.0)
    ((< z 0.0)
     (digamma-negative z))
    ((eq? z 1.0)
     DIGAMMA-1)
    ((<=3D z DIGAMMA-SMALL)
     (/ (- DIGAMMA-1 1.0)
        (+ z (* DIGAMMA-2 z))))
    (else
     (let-values (((z z-shift) (digamma-shift z))) ;; Reduce x < large to
digamma(x+n) where (x+n) >=3D large
       (let ((r (/ z)))                            ;; Use de Moivre's
expansion if argument >=3D large.
         (let ((A (+ z-shift (- (fllog z) (* 0.5 r)))))
           (+ A (digamma-poly (* r r)))))))))


(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
         (=CE=BB: ([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 30)
(plot-error digamma digamma* 0 30)

(plot-error flpsi0 digamma* -8 0)
(plot-error digamma digamma* -8 0)



On Tue, May 7, 2013 at 11:39 AM, Neil Toronto <neil.toronto at gmail.com>wrote=
:

> On 05/06/2013 01:25 PM, Ray Racine wrote:
>
>> Neil,
>>
>> Here in is an implementation of the digamma function.
>>
>
> 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.
>
>
>  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
>> http://code.google.com/p/**pmtk3/source/browse/trunk/**
>> util/digamma.m?r=3D227&spec=3D**svn227<http://code.google.com/p/pmtk3/so=
urce/browse/trunk/util/digamma.m?r=3D227&spec=3Dsvn227>
>> from pmtk2 a probabilistic modeling toolkit for Matlab/Octave, which is
>> MIT LICENSED so I think we are good there.
>>
>
> 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
>
>
>  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.
>>
>
> 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.
>
>
>  Since it works for x =3D 10 to a few decimal places, I'm sure it works
>> everywhere ... perfectly.
>>
>
> :D
>
>
>  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.
>>
>
> 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 <=3D 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 <=3D 0.5 ulps is *correct*. I call a function with =
<=3D
> 2 ulps friggin' fantastic, and one with <=3D 10 accurate.
>
>  I added addition poly coefficients up to 12 but the impl currently only
>> uses 7.
>>
>> impl follows:  [...]
>>
>
> 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
>          (=CE=BB: ([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. <=3D 2 ulps, which is friggin' fantastic). If using 12 terms instea=
d 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 =E2=8A=A5
>
> ____________________
>  Racket Users list:
>  http://lists.racket-lang.org/**users <http://lists.racket-lang.org/users=
>
>

--047d7b3a9684d43f1e04dc27158c
Content-Type: text/html; charset=UTF-8
Content-Transfer-Encoding: quoted-printable

<div dir=3D"ltr"><div><div>Well since I gone this far, I took the unusual s=
tep of implementing the algorithm (more?) correctly this time.<br><br>On th=
e Toronto scale, we are mostly accurate for most of the domain=C2=A0 x &gt;=
 0.=C2=A0 <br>
<br>For x&lt;0 the use of psi(1 - x) - psi (x) =3D pi * cot(pi * x), not so=
 good.<br><br></div>For x &gt; 0, outside the subdomain around the circle x=
=3D1.48 where things are a bit wobbly it is pretty solid.=C2=A0 I don&#39;t=
 know if this is a zero of the polynomial or what the heck is going on ther=
e.<br>
<br></div><div>;; Improved Impl of digamma<br></div><div><div><br>#lang typ=
ed/racket/base<br><br>#|<br>Reference:<br><br>J Bernardo,<br>Psi ( Digamma =
) Function,<br>Algorithm AS 103,<br>Applied Statistics,<br>Volume 25, Numbe=
r 3, pages 315-317, 1976.<br>
<br>From <a href=3D"http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f"=
>http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f</a><br>|#<br><br><b=
r>;; Assumes for z &lt; 0<br>;;=C2=A0 Use the reflection formula (Jeffrey 1=
1.1.6):<br>
;;=C2=A0 digamma(-x) =3D digamma(x+1) + pi*cot(pi*x)<br>;;=C2=A0=C2=A0 y =
=3D digamma(-x+1) + pi*cot(-pi*x);<br>;;=C2=A0 This is related to the ident=
ity<br>;;=C2=A0 digamma(-x) =3D digamma(x+1) - digamma(z) + digamma(1-z)<br=
>;;=C2=A0 where z is the fractional part of x<br>
;;=C2=A0 For example:<br>;;=C2=A0 digamma(-3.1) =3D 1/3.1 + 1/2.1 + 1/1.1 +=
 1/0.1 + digamma(1-0.1)<br>;;=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 =3D digamma(4.1) - digamma(0.=
1) + digamma(1-0.1)<br>;;=C2=A0 Then we use<br>;;=C2=A0 digamma(1-z) - diga=
mma(z) =3D pi*cot(pi*z)<br>
;;=C2=A0 cot(x) =3D tan (pi/2 - x)<br><br>(provide:<br>=C2=A0[digamma (Floa=
t -&gt; Float)])<br><br>(require<br>=C2=A0racket/match<br>=C2=A0racket/pret=
ty<br>=C2=A0(only-in racket/flonum<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0 fltan<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0 fllog)<br>=C2=A0(only-in racket/math<br>
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 pi))<br><br>(define:=
 PI : Float pi)<br>(define: NEG-PI : Float (- PI))<br>(define: PI-OVER-TWO =
: Float (/ PI 2.0))<br>(define: PI-SQUARED-OVER-SIX : Float (/ (* PI PI) 6)=
)<br>(define: EULER-MASCHERONI : Float -0.5772156649015328606065121)<br>
<br>(define: DIGAMMA-1 : Float EULER-MASCHERONI)=C2=A0=C2=A0=C2=A0=C2=A0 ;;=
 digamma(1)<br>(define: DIGAMMA-2 : Float PI-SQUARED-OVER-SIX)<br><br>(defi=
ne: DIGAMMA-COEF-2 : Float (- (/ 1.0 12.0)))<br>(define: DIGAMMA-COEF-3 : F=
loat (/ 1.0 120.0))<br>
(define: DIGAMMA-COEF-4 : Float (- (/ 1.0 252.0)))<br>(define: DIGAMMA-COEF=
-5 : Float (/ 1.0 240.0))<br>(define: DIGAMMA-COEF-6 : Float (- (/ 1.0 132.=
0)))<br>(define: DIGAMMA-COEF-7 : Float (/ 691.0 32760.0))<br>(define: DIGA=
MMA-COEF-8 : Float (- (/ 1.0 12.0)))<br>
<br>(define: DIGAMMA-LARGE : Float 9.5)<br>(define: DIGAMMA-SMALL : Float 0=
.000001)<br><br>;;anxn + =E2=80=A6 + a2x2 + a1x + a0 =3D a0 + x(a1 + x(a2 +=
 x(=E2=80=A6 an) =E2=80=A6) )<br>(: digamma-poly (Float -&gt; Float))<br>(d=
efine (digamma-poly z)<br>
=C2=A0 (* z (+ DIGAMMA-COEF-2 <br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-3<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (* z =
(+ DIGAMMA-COEF-4<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-5<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-6<br>
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-7<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-8))))))))))=
)))))<br><br>(: cotan (Float -&gt; Float))<br>(define (cotan x)<br>=C2=A0 (=
fltan (- PI-OVER-TWO x)))<br>
<br>(: digamma-negative (Float -&gt; Float))<br>(define (digamma-negative z=
)<br>=C2=A0 ;; (assert (&lt; z 0.0))<br>=C2=A0 (+ (digamma (+ (- z) 1.0))<b=
r>=C2=A0=C2=A0=C2=A0=C2=A0 (* PI (cotan (* z NEG-PI)))))<br><br>(: digamma-=
shift (Float -&gt; (Values Float Float)))<br>
(define (digamma-shift z)<br>=C2=A0 (let loop ((z z) (y 0.0))<br>=C2=A0=C2=
=A0=C2=A0 (if (&gt;=3D z DIGAMMA-LARGE)<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0 (values z y)<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (lo=
op (+ z 1.0) (- y (/ z))))))<br><br>(: digamma (Float -&gt; Float))<br>(def=
ine (digamma z)<br>
=C2=A0 (cond <br>=C2=A0=C2=A0=C2=A0 ((eq? z 0.0)<br>=C2=A0=C2=A0=C2=A0=C2=
=A0 -inf.0) <br>=C2=A0=C2=A0=C2=A0 ((&lt; z 0.0)<br>=C2=A0=C2=A0=C2=A0=C2=
=A0 (digamma-negative z))<br>=C2=A0=C2=A0=C2=A0 ((eq? z 1.0)<br>=C2=A0=C2=
=A0=C2=A0=C2=A0 DIGAMMA-1)<br>=C2=A0=C2=A0=C2=A0 ((&lt;=3D z DIGAMMA-SMALL)=
<br>=C2=A0=C2=A0=C2=A0=C2=A0 (/ (- DIGAMMA-1 1.0)<br>=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (+ z (* DIGAMMA-2 z))))<br>
=C2=A0=C2=A0=C2=A0 (else=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 <br>=
=C2=A0=C2=A0=C2=A0=C2=A0 (let-values (((z z-shift) (digamma-shift z))) ;; R=
educe x &lt; large to digamma(x+n) where (x+n) &gt;=3D large=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0 <br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (let ((r (/ z)))=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0 ;; Use de Moivre&#39;s expansion if argument &gt;=3D large.<br>
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (let ((A (+ z-shift (- (fl=
log z) (* 0.5 r)))))<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0 (+ A (digamma-poly (* r r)))))))))<br><br><br>(require plot/typed=
<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 math)<br><br>(: digamm=
a* (Flonum -&gt; Real))<br>(define (digamma* x)<br>
=C2=A0 (parameterize ([bf-precision=C2=A0 53])<br>=C2=A0=C2=A0=C2=A0 (bigfl=
oat-&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>
=C2=A0 (plot (function<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 =
(=CE=BB: ([x : Real])<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0 (let ([x=C2=A0 (fl x)])<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (define err (flulp-error (f-appx x)=
 (f-truth x)))<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0 (if (rational? err) err -1.0)))<br>=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0 x-min x-max)))<br>
<br>(plot-error flpsi0 digamma* 0 30)<br>(plot-error digamma digamma* 0 30)=
<br><br>(plot-error flpsi0 digamma* -8 0)<br>(plot-error digamma digamma* -=
8 0)<br><br></div></div></div><div class=3D"gmail_extra"><br><br><div class=
=3D"gmail_quote">
On Tue, May 7, 2013 at 11:39 AM, Neil Toronto <span dir=3D"ltr">&lt;<a href=
=3D"mailto:neil.toronto at gmail.com" target=3D"_blank">neil.toronto at gmail.com=
</a>&gt;</span> wrote:<br><blockquote class=3D"gmail_quote" style=3D"margin=
:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class=3D"im">On 05/06/2013 01:25 PM, Ray Racine wrote:<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #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=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
Initially I tried<br>
Haskell&#39;s, then a Java impl from Mallet =C2=A0(LDA library) but they di=
dn&#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=3D"http://code.google.com/p/pmtk3/source/browse/trunk/util/digamma.=
m?r=3D227&amp;spec=3Dsvn227" target=3D"_blank">http://code.google.com/p/<u>=
</u>pmtk3/source/browse/trunk/<u></u>util/digamma.m?r=3D227&amp;spec=3D<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 c=
lass=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #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 positiv=
e<br>
x values.<br>
</blockquote>
<br></div>
You can use `bfpsi0&#39; from `math/bigfloat&#39; to get arbitrarily accura=
te error measurements. I used it to produce error plots when implementing `=
flpsi0&#39;. Demo below.<div class=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
Since it works for x =3D 10 to a few decimal places, I&#39;m sure it works<=
br>
everywhere ... perfectly.<br>
</blockquote>
<br></div>
:D<div class=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #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. =C2=A0The 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 func=
tion implementation acceptable when:<br>
<br>
=C2=A01. It gives sensible answers on the largest possible domain. (So<br>
=C2=A0 =C2=A0 Octave+GSL&#39;s digamma doesn&#39;t qualify.)<br>
<br>
=C2=A02. Its error in ulps is &lt;=3D 10 over most of the domain.<br>
<br>
=C2=A03. When error in ulps is higher in some subdomain, it&#39;s because t=
here<br>
=C2=A0 =C2=A0 are no reasonably fast algorithms for computing the function<=
br>
=C2=A0 =C2=A0 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 tru=
e answer and the approximation.<br>

<br>
A function with error &lt;=3D 0.5 ulps is *correct*. I call a function with=
 &lt;=3D 2 ulps friggin&#39; fantastic, and one with &lt;=3D 10 accurate.<b=
r>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex"><div class=3D"im">
I added addition poly coefficients up to 12 but the impl currently only<br>
uses 7.<br>
<br></div>
impl follows: =C2=A0[...]<br>
</blockquote>
<br>
Paste this at the bottom to test it:<br>
<br>
<br>
(require plot/typed<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0math)<br>
<br>
(: digamma* (Flonum -&gt; Real))<br>
(define (digamma* x)<br>
=C2=A0 (parameterize ([bf-precision =C2=A053])<br>
=C2=A0 =C2=A0 (bigfloat-&gt;real (bfpsi0 (bf x)))))<br>
<br>
(: plot-error ((Flonum -&gt; Flonum) (Flonum -&gt; Real) Real Real -&gt; An=
y))<br>
(define (plot-error f-appx f-truth x-min x-max)<br>
=C2=A0 (plot (function<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(=CE=BB: ([x : Real])<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(let ([x =C2=A0(fl x)])<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(define err (flulp-error (f=
-appx x) (f-truth x)))<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(if (rational? err) err -1.=
0)))<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0x-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;=3D 2 ulps, which is friggin&#39; fantastic). If using 1=
2 terms instead of 7 in your polynomial approximation reduced the error for=
 smaller x, it would be worth using that polynomial approximation in the ma=
th library. (I can&#39;t test this because I don&#39;t know the coefficient=
s&#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 slo=
w to evaluate as similar-degree polynomials using Horner&#39;s rule (as in =
your `digamma-poly&#39;).<br>

<br>
Neil =E2=8A=A5<br>
<br>
____________________<br>
=C2=A0Racket Users list:<br>
=C2=A0<a href=3D"http://lists.racket-lang.org/users" target=3D"_blank">http=
://lists.racket-lang.org/<u></u>users</a><br>
</blockquote></div><br></div>

--047d7b3a9684d43f1e04dc27158c--

Posted on the users mailing list.