gammarnd Function

recursive function gammarnd(shape,scale) result(ans)
  real(8) ::  shape,scale,u,w,d,c,x,xsq,g,v,ans
  if (shape <= 0d0) then
     write(*,*) "GAMMARND: Shape parameter must be positive"
  end if
  if (scale <= 0d0) then
     write(*,*) "GAMMARND: Scale parameter must be positive"
  end if
  !    ## Implementation based on "A Simple Method for Generating Gamma Variables"
  !    ## by George Marsaglia and Wai Wan Tsang.  
  !    ## ACM Transactions on Mathematical Software
  !    ## Vol 26, No 3, September 2000, pages 363-372.
  if (shape >= 1d0) then
     d = shape - 1d0/3d0
     c = 1d0/(9d0*d)**0.5d0
     do while (.true.)
        x = normalrnd(0d0,1d0)
        v = 1.d0 + c*x
        do while (v <= 0d0) 
           x = normalrnd(0d0,1d0)
           v = 1d0 + c*x
        end do
        v = v*v*v
        u = mersenne()
        xsq = x*x
        if ((u < 1d0 -.0331d0*xsq*xsq) .or. (log(u) < 0.5d0*xsq + d*(1d0 - v + log(v))) )then
           ans=scale*d*v
           return 
        end if
     end do
  else
     g = gammarnd(shape+1d0,1d0)
     w = mersenne()
     ans=scale*g*(w)**(1d0/shape)
     return 
  end if
end function gammarnd