      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


   
