# Gamma and Polygamma functions

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

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

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

No need to go to the library, Numerical Recipes is all online:

http://www.numerical-recipes.com/nronline_switcher.html

Ben

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

Hello peter, is this what you’re looking for ?

digamma < - function(z)

{

# 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 < - is.na(z) | Re(z) <= 0)) {

ps < - z

ps[omit] < - NA

if(any(!omit))

ps[!omit] < - Recall(z[!omit])

return(ps)

}

if(any(small < - Mod(z) < 5)) {

ps < - z

x < - z[small]

ps[small] < - Recall(x + 5) - 1/x - 1/(x + 1) - 1/(x + 2) - 1/(x + 3) -

1/(x + 4)

if(any(!small))

ps[!small] < - Recall(z[!small])

return(ps)

}

x < - 1/z^2

tail < - ((x * (-1/12 + ((x * (1/120 + ((x * (-1/252 + ((

x * (1/240 + ((x * (-1/132 + ((x * (691/32760 + (

(x * (-1/12 + (3617 * x)/8160)))))))))))))))))))))

log(z) – 1/(2 * z) + tail

}

trigamma < - function(z)

{

# 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 < - is.na(z) | Re(z) <= 0)) {

ps < - z

ps[omit] < - NA

if(any(!omit))

ps[!omit] < - Recall(z[!omit])

return(ps)

}

if(any(small < - Mod(z) < 5)) {

ps < - z

x < - z[small]

ps[small] < - Recall(x + 5) + 1/x^2 + 1/(x + 1)^2 +

1/(x + 2)^2 + 1/(x + 3)^2 + 1/(x + 4)^2

if(any(!small))

ps[!small] < - Recall(z[!small])

return(ps)

}

x < - 1/z^2

tail < - 1 + (x * (1/6 + (x * (-1/30 + (x * (1/42 + (x * (

-1/30 + (x * (5/66 + (x * (-691/2370 + (x * (7/6 -

(3617 * x)/510))))))))))))))

1/(2 * z^2) + tail/z

}

logmdig < - function(z)

{

# log(z) – digamma(z)

# Saves computation of log(z) and avoids subtractive cancellation when z

is large

# GKS 19 Jan 98

#

if(any(omit < - is.na(z) | Re(z) <= 0)) {

ps < - z

ps[omit] < - NA

if(any(!omit))

ps[!omit] < - Recall(z[!omit])

return(ps)

}

if(any(small < - Mod(z) < 5)) {

ps < - z

x < - z[small]

ps[small] < - log(x/(x+5)) + Recall(x + 5) + 1/x + 1/(x + 1) + 1/(x + 2) +

1/(x + 3) + 1/(x + 4)

if(any(!small))

ps[!small] < - Recall(z[!small])

return(ps)

}

x < - 1/z^2

tail < - ((x * (-1/12 + ((x * (1/120 + ((x * (-1/252 + ((

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/

========================

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

> 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.

http://sun.cwru.edu/~jiayang/s.html

Ben

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