function gammam ( x )
       double precision gammam,x
       double precision,external :: gammln
       !x=12.35020017623901
	
       gammam=exp(gammln(x))
      ! print *, "gamma=",gammam
	
       return 
       end