!recursive double precision function Boys(x,n) result(res)
! implicit none
! include 'constants.F'
!
! real, intent(in)     :: x
! integer, intent(in)  :: n
!
! ASSERT (x > 0.)
! if (n == 0) then
!   res = sqrt(pi/(4.*x))*erf(sqrt(x))
! else
!   res = (dble(2*n-1) * Boys(x,(n-1)) - exp(-x) )/(2.*x)
! endif
!
!end function


double precision function binom(k,n)
 implicit none
 integer, intent(in) :: k,n
 double precision, save :: memo(0:20,0:20)
 !DEC$ ATTRIBUTES ALIGN : 32 :: memo
 integer, save :: ifirst
 double precision :: fact
 if (ifirst /= 0) then

  if (n<20) then
    binom = memo(k,n)
    return
  else
    binom=fact(n)/(fact(k)*fact(n-k))
    return
  endif

 else

  ifirst = 1
  integer :: i,j
  do j=0,20
   do i=0,j
    memo(i,j) = fact(j)/(fact(i)*fact(j-i))
   enddo
  enddo

  if (n<20) then
    binom = memo(k,n)
    return
  else
    binom=fact(n)/(fact(k)*fact(n-k))
  endif
 endif
end

double precision function fact2(n)
  implicit none
  integer :: n
  double precision, save :: memo(100)
  integer, save :: memomax = 1 
  
  ASSERT (mod(n,2) /= 0)
  memo(1) = 1.d0
  if (n<2) then
      fact2 = 1.d0
  else if (n<=memomax) then
      fact2 = memo(n)
  else

    integer :: i
    do i=memomax+2,min(n,99),2
      memo(i) = memo(i-2)* float(i)
    enddo
    memomax = min(n,99)
    fact2 = memo(memomax)

    do i=101,n,2
      fact2 *= float(i)
    enddo

  endif

end function


double precision function fact(n)
  implicit none
  integer :: n
  double precision, save :: memo(0:100)
  integer, save :: memomax = 1 
  
  memo(0) = 1.d0
  memo(1) = 1.d0

  if (n<=memomax) then
    fact = memo(n)
  else
    integer :: i
    do i=memomax+1,min(n,100)
      memo(i) = memo(i-1)*float(i)
    enddo
    memomax = min(n,100)
    fact = memo(memomax)
    do i=101,n
      fact *= float(i)
    enddo
  endif
end function


double precision function rintgauss(n)
  implicit none
  include 'constants.F'

  integer :: n

  rintgauss = dsqrt_pi
  if ( n == 0 ) then
    return
  else if ( n == 1 ) then
    rintgauss = 0.
  else if ( mod(n,2) == 1) then
    rintgauss = 0.
  else
    double precision :: fact2
    rintgauss = rintgauss/(2.**(n/2))
    rintgauss = rintgauss * fact2(n-1)
  endif
end function

double precision function goverlap(gamA,gamB,nA)
  implicit none

  real    :: gamA, gamB
  integer :: nA(3)

  double precision :: gamtot
  gamtot = gamA+gamB

  goverlap=1.0

  integer :: l
  double precision :: rintgauss
  do l=1,3
    goverlap = goverlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
  enddo

end function