      function beta(a,b) ! bernd berg apr 02 2000.
c incomplete beta function
      include '../../libs/fortran/implicit.sta'
      beta=exp(gamma_ln(a)+gamma_ln(b)-gamma_ln(a+b))
      return
      end
