Forums > Dev

Gamma and Polygamma functions

February 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


February 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


February 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

http://www.bloghotel.org/LIEU/


February 14, 2006 | 3:38 pm

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

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

Ben


February 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


February 15, 2006 | 2:58 pm

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/
========================


February 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


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

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

Ben


February 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


Viewing 9 posts - 1 through 9 (of 9 total)