Gamma and Polygamma functions


    Feb 14 2006 | 9:57 am
    Dear Collected Wisdom,
    Does anyone here know of any available Math libraries that implement
    the polygamma functions (digamma, trigamma, etc., usually notated as
    psi)? I've not had success with/Google (tried both web and newsgroups).
    Also, I'm not clear whether CodeWarrior implements the ordinary
    gamma/log-gamma functions on Windows, so I may need some code for that
    as well. I think Numerical Recipes has sample code, I'll try to check
    that at the library.
    Thanks for ideas/pointers/etc.
    Best -- Peter
    -------------- http://www.bek.no/~pcastine/Litter/ --------------
    Peter Castine | ^
    | Litter Power & Litter Bundle for Jitter
    pcastine@gmx.net |
    pcastine@bek.no | iCE: Sequencing, Recording, and Interface Building
    4-15@kagi.com | for Max/MSP
    | Extremely cool
    | http://www.dspaudio.com
    | http://www.dspaudio.com/software/software.html

    • Feb 14 2006 | 10:05 am
      hi peter
      i remember that xenakis ultra core book, formalized music
      has some pointers related to this. but you probably know that already
      hope i'm not confusing what you are asking, and hope this helps
      a
    • Feb 14 2006 | 10:15 am
      I have a great book of that kind of things in C.
      Not here but in Lyon where i have to go back tomorrow.
      I'll have a look at it to see if i can find something that would fit your
      need.
      Talk to you soon then :-)
      Sebowski
    • Feb 14 2006 | 3:38 pm
      No need to go to the library, Numerical Recipes is all online:
      Ben
    • Feb 15 2006 | 10:06 am
      Thanks to all for the pointers to NR online.
      Funnily enough, the code for lgamma() made a lot more sense to me when
      I was holding the dead-tree copy in my hands. OTOH, I'm glad enough to
      be able to copy&paste the code, particularly the float constants. And
      it seems strange that the Mac math libraries provide lgamma() while
      Windows doesn't, given that all my EE friends live in Windowsland.
      <>.
      Anyway, I'm still in the market for polygamma code. I'll keep on
      looking.
      -- Peter
      -------------- http://www.bek.no/~pcastine/Litter/ --------------
      Peter Castine | ^
      | Litter Power & Litter Bundle for Jitter
      pcastine@gmx.net |
      pcastine@bek.no | iCE: Sequencing, Recording, and Interface Building
      4-15@kagi.com | for Max/MSP
      | Extremely cool
      | http://www.dspaudio.com
      | http://www.dspaudio.com/software/software.html
    • Feb 15 2006 | 2:58 pm
      Hello peter, is this what you're looking for ?
      digamma {
      # Digamma function: first derivative of log(gamma(x)).
      # Original function by Bill Venables.
      # Modified by GKS to accepts NAs, 19 Jan 1998.
      #
      if(any(omit ps ps[omit] if(any(!omit))
      ps[!omit] return(ps)
      }
      if(any(small ps x ps[small] 1/(x + 4)
      if(any(!small))
      ps[!small] return(ps)
      }
      x tail x * (1/240 + ((x * (-1/132 + ((x * (691/32760 + (
      (x * (-1/12 + (3617 * x)/8160)))))))))))))))))))))
      log(z) - 1/(2 * z) + tail
      }
      trigamma {
      # Trigamma function: second derivative of log(gamma(x)).
      # Original function by Bill Venables.
      # Modified by GKS to accepts NAs, 19 Jan 1998.
      #
      if(any(omit ps ps[omit] if(any(!omit))
      ps[!omit] return(ps)
      }
      if(any(small ps x ps[small] 1/(x + 2)^2 + 1/(x + 3)^2 + 1/(x + 4)^2
      if(any(!small))
      ps[!small] return(ps)
      }
      x tail -1/30 + (x * (5/66 + (x * (-691/2370 + (x * (7/6 -
      (3617 * x)/510))))))))))))))
      1/(2 * z^2) + tail/z
      }
      logmdig {
      # log(z) - digamma(z)
      # Saves computation of log(z) and avoids subtractive cancellation when z
      is large
      # GKS 19 Jan 98
      #
      if(any(omit ps ps[omit] if(any(!omit))
      ps[!omit] return(ps)
      }
      if(any(small ps x ps[small] 1/(x + 3) + 1/(x + 4)
      if(any(!small))
      ps[!small] return(ps)
      }
      x tail x * (1/240 + ((x * (-1/132 + ((x * (691/32760 + (
      (x * (-1/12 + (3617 * x)/8160)))))))))))))))))))))
      1/(2 * z) - tail
      }
      http://www.statsci.org/s/polygamm.s
      ========================
      Seb.Tworowski
      http://www.bloghotel.org/LIEU/
      ========================
    • Feb 16 2006 | 12:46 am
      On around Feb 15, 2006, at 15:58, tworowski.sebastien@laposte.net said
      something like:
      > Hello peter, is this what you're looking for ?
      Yes indeed.
      BTW, what language is it in? I can figure out most of the code, but one
      or two things are a little strange.
      Thanks for sharing this!
      Best -- Peter
      >
      -------------- http://www.bek.no/~pcastine/Litter/ --------------
      Peter Castine | ^
      | Litter Power & Litter Bundle for Jitter
      pcastine@gmx.net |
      pcastine@bek.no | iCE: Sequencing, Recording, and Interface Building
      4-15@kagi.com | for Max/MSP
      | Extremely cool
      | http://www.dspaudio.com
      | http://www.dspaudio.com/software/software.html
    • Feb 16 2006 | 5:35 am
      > BTW, what language is it in?
      It's s+, a programming language for data analysis and statistics. I
      used it once or twice in my last job.
      Ben
    • Feb 16 2006 | 9:22 am
      On around Feb 16, 2006, at 6:35, Ben Nevile said something like:
      > It's s+, a programming language for data analysis and statistics.
      Cool.
      > I used it once or twice in my last job.
      Which means that you're someone who would actually understand the math
      inside the Litter Power RNGs. Also cool.
      I managed to work out Digamma for ints and half-ints, but the entire
      real domain needs incantations I don't know. Thanks 2^20.
      -- Peter
      >
      -------------- http://www.bek.no/~pcastine/Litter/ --------------
      Peter Castine | ^
      | Litter Power & Litter Bundle for Jitter
      pcastine@gmx.net |
      pcastine@bek.no | iCE: Sequencing, Recording, and Interface Building
      4-15@kagi.com | for Max/MSP
      | Extremely cool
      | http://www.dspaudio.com
      | http://www.dspaudio.com/software/software.html