program hw02_01 implicit none * Program to generate two tables of the gamma function accurate to * at least significant digits. * Methods are based on: * * M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with * Formulas, Graphs and Mathematical Tables, Wiley, New York, 1970. * (Gamma finctons: page 255) * Web access: * http://www.convertit.com/Go/Convertit/Reference/AMS55.ASP * http://www.convertit.com/Go/Convertit/Reference/AMS55.ASP?Res=150&Page=255 * * Basic methosd here is explicit formualas of Gamma(n) = (n-1)!, * gamma(n+1/3, 1/2 and 2/3). For the -1 to +1 ranges, Euler's series is used * real*8 gammai, gammathird, gammatwothird, gammahalf real*8 z, gammaeul,gammainf, gammaser real*8 eps, gamma_1, gamma_2, gamma_3 character*80 line integer*4 i *** Start program write(*,100) 100 format(/,'Table of gamma functions of positive n integers',/, . '-----------------------------------------------',/, . ' n gamma(n) gamma(n+1/3) gamma(n+1/2)', . ' gamma(n+2/3)') do i = 1,10 write(*,120) i, gammai(i), gammathird(i), gammahalf(i), . gammatwothird(i) 120 format(I4,2x,F12.0,1x,F15.6,2x,F15.6,2x,F15.6) end do eps = 1.d-6 write(*,200) eps 200 format(/,'Table of gamma functions for non-integer arguments',/, . '--------------------------------------------------',/, . 'Three methods are used here:',/, . 'Eulers formula GammaEul',/, . 'Eulers infinite product GammaInf',/, . 'Asymptotic series expansion GammaSer',/, . 'Tolerance on calculation is ',E12.3,/, . ' z GammaEul(z) GammaInf(z) ', . ' GammaSer(z)') do i = -10,10 z = i/10.d0 gamma_1 = gammaeul(z,eps) gamma_2 = gammainf(z,eps) gamma_3 = gammaser(z,eps) *** Formthe output line if( gamma_1.eq. 0 ) then write(line(1:),210) char(177) 210 format(6x,a1,'Infinity' ) else write(line(1:),215) gamma_1 215 format(1x,F14.6) endif if( gamma_2.eq. 0 ) then write(line(21:),210) char(177) else write(line(21:),215) gamma_2 endif if( gamma_3.eq. 0 ) then write(line(41:),210) char(177) else write(line(41:),215) gamma_3 endif write(*,220) z, line(1:61) 220 format(F8.2,2x,a60) enddo end real*8 function gammai(n) implicit none * Function to return gamma(n) when n is positive interger * gamma(n) = (n-1)! * PASSED VARIABLE integer*4 n * LOCAL VARIABLES integer*4 i ! Loop counter while computing factorial * Set initial value and then loop to (n-1) forming product gammai = 1.d0 if( n.lt.1 ) then write(*,100) n 100 format('Invalid argument in gammai call: ', . 'Argument must be > 1. Passed value was ',i6) stop 'Invalid argument to gammai function' endif do i = 2,n-1 gammai = gammai*i end do ***** End return end real*8 function gammathird(n) implicit none * Function to return gamma(n+1/3) when n is positive interger * gamma(n+1/3) = (1*4*7*...*(3n-2)*gamma(1/3)/3**n * PASSED VARIABLE integer*4 n * LOCAL VARIABLES integer*4 i ! Loop counter while computing factorial * Set initial value and then loop to (n-1) forming product gammathird = 2.6789385347077479d0 if( n.lt.0 ) then write(*,100) n 100 format('Invalid argument in gammathird call: ', . 'Argument must be >=0. Passed value was ',i6) stop 'Invalid argument to gammathird function' endif do i = 4,(3*n-2),3 gammathird = gammathird*i end do gammathird = gammathird /(3.d0**n) ***** End return end real*8 function gammatwothird(n) implicit none * Function to return gamma(n+2/3) when n is positive interger * gamma(n+1/3) = (2*5*8*...*(3n-1)*gamma(2/3)/3**n * PASSED VARIABLE integer*4 n * LOCAL VARIABLES integer*4 i ! Loop counter while computing factorial * Set initial value and then loop to (n-1) forming product gammatwothird = 1.3541179394264005d0 if( n.lt.0 ) then write(*,100) n 100 format('Invalid argument in gammatwothird call: ', . 'Argument must be >=0. Passed value was ',i6) stop 'Invalid argument to gammatwothird function' endif do i = 2,(3*n-1),3 gammatwothird = gammatwothird*i end do gammatwothird = gammatwothird /(3.d0**n) ***** End return end real*8 function gammahalf(n) implicit none * Function to return gamma(n+1/2) when n is positive interger * gamma(n+1/2) = (1*3*5*...*(2n-1)*gamma(1/2)/2**n * PASSED VARIABLE integer*4 n * LOCAL VARIABLES integer*4 i ! Loop counter while computing factorial * Set initial value and then loop to (n-1) forming product gammahalf = 1.7724538509055161d0 if( n.lt.0 ) then write(*,100) n 100 format('Invalid argument in gammahalf call: ', . 'Argument must be >= 0. Passed value was ',i6) stop 'Invalid argument to gammahalf function' endif do i = 3,(2*n-1),2 gammahalf = gammahalf*i end do gammahalf = gammahalf /(2.d0**n) ***** End return end real*8 function gammaeul(z,eps) implicit none * Function to return gamma(z) where z is any non-negative integer values * The expansion here is based on Euler's formulas: * gamma(z) = limit(n -> infinity) n! n^z/(z*(z+1)*(z_2)*...(z+n) * for all z except 0 and negative integers. The restriction can be clearly * seen because for these values of z we would divide by zero and gamma at * these valuse is infinite * * To determine the size of n needed the above is re-arranged so that each * term can be added i.e., * gamma_n = gamma_(n-1) *n/(z+n) * n**z/(n-1)**z * by looking at second part we can see if enough terms have been included * real*8 z ! argument ., eps ! Tolerance on accuracy needed real*8 gamma_n, gamma_prev ! Iterated values of gamma, gamma_n and gamma_(n-1) real*8 err ! Error in current iteration real*8 npowz,nm1powz ! n and n-1 to the power of z. Save value so that we do not need to ! re-compute integer*4 n ! Counter ., max_n ! Max n allowed (to stop iteration) max_n = 10000000 if( z.le.0 .and. abs(z-nint(z)).lt.eps ) then c write(*,*) 'z too close to 0,-1,-2, ... ' gammaeul = 0.d0 return endif n = 1 gamma_prev = 1.d0/(z*(z+1)) err = 1.d0 nm1powz = 1.0d0 ! (n-1)^z for n = 1 do while ( n.lt. max_n .and. err.gt.eps/1.d6) n = n + 1 npowz = (dble(n)**z) gamma_n = gamma_prev * (dble(n)/(z+n)) * . (npowz/nm1powz) nm1powz = npowz ! Save value for next iteration err = abs(gamma_n - gamma_prev) gamma_prev = gamma_n end do c write(*,*) 'Gamma ', gamma_n, n, err gammaeul = gamma_n return end real*8 function gammainf(z,eps) implicit none * Function to return gamma(z) where z is any non-negative integer values * The expansion here is based on Euler's Infinit product formula of 1/gamma: * 1/gamma(z) = z*exp(gammacon*z)*product_n[(1+z/n)exp(-z/n)] * for all z. When 1/gamma is zero the gamma is infinite. Gammma_con is the * constant gamma = 0.5772 15664 .. * * To determine the size of n needed is determined by looking at the change * at step. The series is very slow to converge. * real*8 z ! argument ., eps ! Tolerance on accuracy needed real*8 gammacon, gamma_prev ! Iterated values of gamma, gamma_n and gamma_(n-1) real*8 gammainv ! Inverse of gamma real*8 factor real*8 err ! Error in current iteration integer*4 n ! Counter ., max_n ! Max n allowed (to stop iteration) max_n = 10000000 if( z.le.0 .and. abs(z-nint(z)).lt.eps ) then c write(*,*) 'z too close to 0,-1,-2, ... ' gammainf = 0.d0 return endif n = 0 gammacon = 0.5772 15663 90153 218606d0 gamma_prev = z*exp(gammacon*z) err = 1 do while ( n.lt.max_n .and. err.gt.eps/1.d6 ) n = n + 1 gammainv = gamma_prev*(1+z/n)*exp(-z/n) err = abs(gammainv-gamma_prev) gamma_prev = gammainv ! print *,'z, n ',z,n, err, gamma_prev, gammainv, 1.d0/gammainv end do c write(*,*) 'Gamma Inf', gammainv, 1.d0/gammainv, n, err gammainf = 1.d0/gammainv return end real*8 function gammaser(z,eps) implicit none * Function to return gamma(z) where z is any non-negative integer values. * This is a series expansion developed in the 1930s for 1/gamma(z) * The coefficients c k are from H. T. Davis, Tables of * higher mathematical functions, 2 vols., Principia Press, * Bloomington, Ind., 1933, 1935 (with permission) ; with * corrections due to H. E. Salaer. * Values from * http://www.convertit.com/Go/Convertit/Reference/AMS55.ASP?Res=150&Page=256 * real*8 z ! argument ., eps ! Tolerance on accuracy needed real*8 coeff(26) ! 26 coefficients in series expansion ., gamma_inv integer*4 n data coeff / . 1.00000 00000 000000d0, . 0.57721 56649 015329d0, . -0.65587 80715 202538d0, . -0.04200 26350 340952d0, . 0.16653 86113 822915d0, . -0.04219 77345 555443d0, . -0.00962 19715 278770d0, . 0.00721 89432 466630d0, . -0.00116 51675 918591d0, . -0.00021 52416 741149d0, . 0.00012 80502 823882d0, . -0.00002 01348 547807d0, . -0.00000 12504 934821d0, . 0.00000 11330 272320d0, . -0.00000 02056 338417d0, . 0.00000 00061 160950d0, . 0.00000 00050 020075d0, . -0.00000 00011 812746d0, . 0.00000 00001 043427d0, . 0.00000 00000 077823d0, . -0.00000 00000 036868d0, . 0.00000 00000 005100d0, . -0.00000 00000 000206d0, . -0.00000 00000 000054d0, . 0.00000 00000 000014d0, . 0.00000 00000 000001d0 / * Check for non-valid of z (inverse would be zero and infinite in the * value. We return zero if( z.le.0 .and. abs(z-nint(z)).lt.eps ) then c write(*,*) 'z too close to 0,-1,-2, ... ' gammaser = 0.d0 return endif * Sum up the 26 terms in the series. gamma_inv = 0.d0 do n = 1,26 gamma_inv = gamma_inv + coeff(n)*z**n end do gammaser = 1.d0/gamma_inv return end