read(*,*) a,b
      n=1000
      h=(b-a)/n
      ip=0
      do i=0,n
         x=a+h*i
         if (x.eq.int(x)) then
            write(*,*)
         else
            gam=xgamma(x)
            write(*,*) x, gam
         end if
      end do
C
      end
C
C------------------------------------------------------
C
      function xgamma(xx)
C
C x! on [0,1], Hastings approx.
C
      x01(x)=
     *1+x*(- 0.577191652 + x*(0.988205891 +
     *  x*(- 0.897056934 + x*(0.918206857 +
     *  x*(- 0.756704078 + x*(0.482199394 +
     *  x*(- 0.193527818 + x*0.035868343)))))))
C
C x! from Stirling approximation
C
      sx(x)=
     *  exp(1/(156.*x**13) 
     *   - 691/(360360.*x**11) 
     *   + 1/(1188.*x**9) 
     *   - 1/(1680.*x**7)
     *   + 1/(1260.*x**5) 
     *   - 1/(360.*x**3) 
     *   + 1/(12.*x)
     *   - x + Log(2*Pi)/2. + (0.5 + x)*Log(x))
      Pi=4*atan(1.0)
C
C------------------------------------------------------
C
C     ..wyzej podane formuly licza z! a nie Gamma(z)=(z-1)!..
      xsig=sign(1.,xx-1)
      x=abs(xx-1)
C
C     ..x = -integer, the pole..
      if( xx.eq.-int(abs(xx)) ) stop 'in xgamma: x = -integer'
C
      if( x.gt.6 ) go to 2
C
C     ..use of Hastings approx. on the [0,1] interval..
C
C     ..x > 1: x!=x*(x-1)!; reduction to [0,1]
C              and reflection formula for negative xx..
      fac=1
    1 if( x.lt.1 ) then
         xgamma=fac*x01(x)
         if( xsig.eq.-1. ) xgamma=Pi*(xx-1)/sin(Pi*(xx-1))/xgamma
         return
      else
         fac=fac*x
         x=x-1
         go to 1
      end if
C
C     ..Stirling approximation..
C
    2 xgamma=sx(x)
      if( xsig.eq.-1) xgamma=Pi*(xx-1)/sin(Pi*(xx-1))/xgamma
      return
C
      end