!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 fact2(n)
  implicit none
  integer :: n
  double precision, save :: memo(1:100)
  integer, save :: memomax = 1 
  
  ASSERT (mod(n,2) /= 0)
  if (n<=memomax) then
    if (n<3) then
      fact2 = 1.d0
    else
      fact2 = memo(n)
    endif
    return
  endif

  integer :: i
  memo(1) = 1.d0
  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 = fact2*float(i)
  enddo

end function


double precision function fact(n)
  implicit none
  integer :: n
  double precision, save :: memo(1:100)
  integer, save :: memomax = 1 
  
  if (n<=memomax) then
    if (n<2) then
      fact = 1.d0
    else
      fact = memo(n)
    endif
    return
  endif

  integer :: i
  memo(1) = 1.d0
  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 = fact*float(i)
  enddo
end function


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

  integer :: n

  rintgauss = sqrt(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