subroutine gaussian_product(a,xa,b,xb,k,p,xp)
 implicit none
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}

 double precision, intent(in) :: a,b   ! Exponents
 double precision, intent(in) :: xa,xb ! Centers
 double precision, intent(out) :: p    ! New exponent
 double precision, intent(out) :: xp   ! New center
 double precision, intent(out) :: k  ! Constant

 double precision:: p_inv

 p = a+b
 xp = (a*xa+b*xb)

 p_inv = 1./p

 xp = xp*p_inv
 k = dexp(-a*b*p_inv*(xa-xb)**2)

end subroutine

double precision function primitive_overlap(a,xa,i,b,xb,j)
 implicit none
 include 'constants.F'
 double precision, intent(in)    :: a,b   ! Exponents
 double precision, intent(in)    :: xa,xb ! Centers
 integer, intent(in) :: i,j   ! Powers of xa and xb

 integer :: ii, jj, kk, ll
 double precision:: xp
 double precision:: p
 double precision :: S(0:i,0:j)
 double precision :: inv_2p, di(i), dj(j)

 call gaussian_product(a,xa,b,xb,S(0,0),p,xp)
 S(0,0) = S(0,0) * dsqrt(pi/p)

 if (i>0) then
   S(1,0) = (xp-xa) * S(0,0) ! TODO verifier signe
 endif
 if (j>0) then
   S(0,1) = (xp-xb) * S(0,0) ! TODO verifier signe
 endif

 inv_2p = 1./(2.*p)

 do ii=1,max(i,j)
  di(ii) = inv_2p * dble(ii)
 enddo

 if (i>1) then
  do ii=1,i-1
   S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
  enddo
 endif

 if (j>1) then
  do jj=1,j-1
   S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
  enddo
 endif

 do jj=1,j
  do ii=1,i
   S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
  enddo
 enddo

 primitive_overlap = S(i,j)

end function