| |
 |
|
|
Science Forum Index » Math - Numerical Analysis Forum » increasing width
Page 1 of 1
|
| Author |
Message |
| Lane Straatman |
 Posted: Fri Mar 16, 2007 8:13 pm |
|
|
|
Guest
|
program gamma2
IMPLICIT NONE
integer, parameter :: DP = SELECTED_REAL_KIND(13,37)
real :: x
real(KIND=DP) :: test
x= 5.7
test = exp(gammln(x))
write (*, '(a, f12. ') 'test = ', test
write (*, '(a)') 'test = 72.527634520222929046'
contains
FUNCTION gammln(xx)
!from Numerical Recipes 6.1
IMPLICIT NONE
REAL gammln, xx
INTEGER j
real(KIND=DP) :: ser, stp, tmp, x, y, cof(6)
! SAVE cof, stp
DATA cof, stp/76.18009172947146d0, -86.5053203941677d0, &
& 24.01409824083091d0, -1.231739572450155d0, &
& .1208650973866179d-2, -.5395239384953d-5,&
& 2.5066282746310005d0/
x=xx
y=x
tmp=x+5.5d0
tmp=(x+0.5d0)*log(tmp)-tmp
ser=1.000000000190015d0
do j=1,6
y=y+1.d0
ser=ser+cof(j)/y
end do
gammln=tmp+log(stp*ser/x)
return
END FUNCTION
end program gamma2
I've been told, in my first thread in this forum, to avoid factorial. To my
thinking, that means that a gamma function is something I must then put in
my grabbag of functions. Right now I get 6 sigfigs with this log gamma
function. I'd like to stretch that with as many different tricks as I can
summon. Suggestions?
--
LS |
|
|
| Back to top |
|
| user923005 |
| Posted: Sat Mar 17, 2007 1:19 am |
|
|
|
Guest
|
On Mar 16, 6:13 pm, "Lane Straatman" <inva...@invalid.net> wrote:
Quote: program gamma2
IMPLICIT NONE
integer, parameter :: DP = SELECTED_REAL_KIND(13,37)
real :: x
real(KIND=DP) :: test
x= 5.7
test = exp(gammln(x))
write (*, '(a, f12. ') 'test = ', test
write (*, '(a)') 'test = 72.527634520222929046'
contains
FUNCTION gammln(xx)
!from Numerical Recipes 6.1
IMPLICIT NONE
REAL gammln, xx
INTEGER j
real(KIND=DP) :: ser, stp, tmp, x, y, cof(6)
! SAVE cof, stp
DATA cof, stp/76.18009172947146d0, -86.5053203941677d0, &
& 24.01409824083091d0, -1.231739572450155d0, &
& .1208650973866179d-2, -.5395239384953d-5,&
& 2.5066282746310005d0/
x=xx
y=x
tmp=x+5.5d0
tmp=(x+0.5d0)*log(tmp)-tmp
ser=1.000000000190015d0
do j=1,6
y=y+1.d0
ser=ser+cof(j)/y
end do
gammln=tmp+log(stp*ser/x)
return
END FUNCTION
end program gamma2
I've been told, in my first thread in this forum, to avoid factorial. To my
thinking, that means that a gamma function is something I must then put in
my grabbag of functions. Right now I get 6 sigfigs with this log gamma
function. I'd like to stretch that with as many different tricks as I can
summon. Suggestions?
--
LS
Whoever it was who told you to avoid factorial was out of their gord.
With 128 bit floating point, there are less than 2000 entries in the
table of answers.
So compute them before-hand, poke them in a table, and look them up.
In the case of 8 byte IEEE floating point, there are a paultry 170 of
them that fit into the table before the exponent overflows.
static const double factorials[] = {
1.0,
1.0,
2.0,
6.0,
2.4e+001,
1.2e+002,
7.2e+002,
5.04e+003,
4.032e+004,
3.6288e+005,
3.6288e+006,
3.99168e+007,
4.790016e+008,
6.2270208e+009,
8.71782912e+010,
1.307674368e+012,
2.0922789888e+013,
3.55687428096e+014,
6.402373705728e+015,
1.21645100408832e+017,
2.43290200817664e+018,
5.109094217170944e+019,
1.1240007277776077e+021,
2.5852016738884978e+022,
6.2044840173323941e+023,
1.5511210043330986e+025,
4.0329146112660565e+026,
1.0888869450418352e+028,
3.0488834461171384e+029,
8.8417619937397008e+030,
2.6525285981219103e+032,
8.2228386541779224e+033,
2.6313083693369352e+035,
8.6833176188118859e+036,
2.9523279903960412e+038,
1.0333147966386144e+040,
3.7199332678990118e+041,
1.3763753091226343e+043,
5.2302261746660104e+044,
2.0397882081197442e+046,
8.1591528324789768e+047,
3.3452526613163803e+049,
1.4050061177528798e+051,
6.0415263063373834e+052,
2.6582715747884485e+054,
1.1962222086548019e+056,
5.5026221598120885e+057,
2.5862324151116818e+059,
1.2413915592536073e+061,
6.0828186403426752e+062,
3.0414093201713376e+064,
1.5511187532873822e+066,
8.0658175170943877e+067,
4.2748832840600255e+069,
2.3084369733924138e+071,
1.2696403353658276e+073,
7.1099858780486348e+074,
4.0526919504877221e+076,
2.3505613312828789e+078,
1.3868311854568986e+080,
8.3209871127413916e+081,
5.0758021387722484e+083,
3.1469973260387939e+085,
1.9826083154044401e+087,
1.2688693218588417e+089,
8.2476505920824715e+090,
5.4434493907744307e+092,
3.6471110918188683e+094,
2.4800355424368305e+096,
1.711224524281413e+098,
1.197857166996989e+100,
8.5047858856786218e+101,
6.1234458376886077e+103,
4.4701154615126834e+105,
3.3078854415193856e+107,
2.4809140811395391e+109,
1.8854947016660498e+111,
1.4518309202828584e+113,
1.1324281178206295e+115,
8.9461821307829729e+116,
7.1569457046263779e+118,
5.7971260207473655e+120,
4.7536433370128398e+122,
3.9455239697206569e+124,
3.314240134565352e+126,
2.8171041143805494e+128,
2.4227095383672724e+130,
2.1077572983795269e+132,
1.8548264225739836e+134,
1.6507955160908452e+136,
1.4857159644817607e+138,
1.3520015276784023e+140,
1.24384140546413e+142,
1.1567725070816409e+144,
1.0873661566567424e+146,
1.0329978488239052e+148,
9.916779348709491e+149,
9.6192759682482062e+151,
9.426890448883242e+153,
9.3326215443944096e+155,
9.3326215443944102e+157,
9.4259477598383536e+159,
9.6144667150351211e+161,
9.9029007164861754e+163,
1.0299016745145622e+166,
1.0813967582402903e+168,
1.1462805637347078e+170,
1.2265202031961373e+172,
1.3246418194518284e+174,
1.4438595832024928e+176,
1.5882455415227421e+178,
1.7629525510902437e+180,
1.9745068572210728e+182,
2.2311927486598123e+184,
2.5435597334721862e+186,
2.9250936934930141e+188,
3.3931086844518965e+190,
3.969937160808719e+192,
4.6845258497542883e+194,
5.5745857612076033e+196,
6.6895029134491239e+198,
8.09429852527344e+200,
9.8750442008335976e+202,
1.2146304367025325e+205,
1.5061417415111404e+207,
1.8826771768889254e+209,
2.3721732428800459e+211,
3.0126600184576582e+213,
3.8562048236258025e+215,
4.9745042224772855e+217,
6.4668554892204716e+219,
8.4715806908788174e+221,
1.1182486511960039e+224,
1.4872707060906852e+226,
1.9929427461615181e+228,
2.6904727073180495e+230,
3.6590428819525472e+232,
5.0128887482749898e+234,
6.9177864726194859e+236,
9.6157231969410859e+238,
1.346201247571752e+241,
1.8981437590761701e+243,
2.6953641378881614e+245,
3.8543707171800706e+247,
5.5502938327393013e+249,
8.0479260574719866e+251,
1.1749972043909099e+254,
1.7272458904546376e+256,
2.5563239178728637e+258,
3.8089226376305671e+260,
5.7133839564458505e+262,
8.6272097742332346e+264,
1.3113358856834518e+267,
2.0063439050956811e+269,
3.0897696138473489e+271,
4.7891429014633912e+273,
7.4710629262828905e+275,
1.1729568794264138e+278,
1.8532718694937338e+280,
2.9467022724950369e+282,
4.714723635992059e+284,
7.5907050539472148e+286,
1.2296942187394488e+289,
2.0044015765453015e+291,
3.2872185855342945e+293,
5.423910666131586e+295,
9.0036917057784329e+297,
1.5036165148649983e+300,
2.5260757449731969e+302,
4.2690680090047027e+304,
7.257415615307994e+306
};
P.S.
There are tons of good gamma() functions all over the net. It makes
no sense to write your own.
Cephes is a good spot to check for a C version.
Here is a Fortran version:
http://ftp.cac.psu.edu/pub/ger/fortran/hdk/gamma.for
Here is another (for log(gamma(x))):
http://lib.stat.cmu.edu/apstat/245
I did not bother to test them (I prefer C or C++ most of the time). |
|
|
| Back to top |
|
| |
|
Page 1 of 1
All times are GMT - 5 Hours
The time now is Sun Sep 07, 2008 10:49 am
|
|