/************************************************

  The function FCint() calculates the overlap integral

  <n1 n2 n3 n4 ...nN | n1' n2' n3' .... nN'>

  that is recursively defined as
  <n1 n2 n3 n4 .ni+1 ..nN | n1' n2' n3' .... nN'> =
     sum_{j=1}^N' 2*sqrt[{nj'}/{ni+1}] RJdagg[i][j] 
         <n1 n2 n3 n4 .ni ..nN | n1' n2' n3' .nj'-1... nN'>
     + sum_{k=1}^N sqrt[{nk}/{ni+1}] twoR1[i][k] 
         <n1 n2 n3 n4 .nk-1 ..nN | n1' n2' n3' .. nN'>
     - sqrt[2/{ni+1}] RJD[i]
         <n1 n2 n3 n4 .ni ..nN | n1' n2' n3' .. nN'>
  to reduce one of the indices of the vector of the n[], and
  recursively defined as
  <n1 n2 n3 n4 ...nN | n1' n2' n3' ..nj'+1.. nN'> =
     sum_{i=1}^N 2*sqrt[{ni}/{nj'+1}] JR[j][i] 
         <n1 n2 n3 n4 .ni-1 ..nN | n1' n2' n3' .. nN'>
     + sum_{k=1}^N' sqrt[{nk'}/{nj'+1}] twoP1[j][k] 
         <n1 n2 n3 n4 ..nN | n1' n2' n3' .nk'-1... nN'>
     + sqrt[2/{nj'+1}] onePD[j]
         <n1 n2 n3 n4 .ni ..nN | n1' n2' n3'.. nj'.. nN'>
  to reduce one of the indices of the vector of the n'[].
  
  The value  of a <n1 n2 n3 n4 ...nN | n1' n2' n3' .... nN'>
  with any of the N numbers n[] or any of the N' number n'[] less than zero
  is defined to be zero.

  Care has been take that the interface remains valid if called
  from Fortran and linked with Fortran programs:
  i) all arguments are pointers
  ii) there are no arguments pointing to string or boolean variables
  iii) the program can be compiled with the macro FCINT_USE_TRANSPOSED_MATRICES
       defined which internally transposes the matrix representations.

  This implements eq 8 and 9 of the article by K K Lehman
  (http://www.princeton.edu/~lehmann) entitled "Recursion Relationships for
  Harmonic Franck-Condon Factors in Polyatomic Molecules" in the version
  of 03-09-1996, which is a slight generalization of eq 11 and 12 of
  M Klessinger, "Calculation of the Vibronic Fine Structure in Electronic
  Spectra at Higher Temperatures. 1. Benzene and Pyrazine", J. Phys. Chem A 102
  (36) 7157-67 (1998)

  The C++ interface definition is

  extern double FCint(const int *n, const int &N, const int *nprime, const int &Nprime, const double *RJdagg, const double *twoR1,
		const double *RJD, const double *JR, const double *twoP1, const double *onePD,const double &OO, const int &init) ;

  Tested with Forte Developer 6 update 2 C++ compiler V 5.1. To compile in this case, use
	CC -c FCint.c
  Richard Mathar, http://www.strw.leidenuniv.nl/~mathar, 21 Aug 2001
***********************************************************/

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <iostream.h>
// <map>, <vector> as contained in the STL
#include <map>
#include <vector>

using namespace std ;

// If called from Fortran, a matrix transposition usually is needed. Instead
// of doing this once before calling FCint() in the Fortran application program, the same results are achieved by swapping
// matrix indices on the fly within FCint(), which is triggered by compilation with FCINT_USE_TRANSPOSED_MATRICES defined.
#define FCINT_USE_TRANSPOSED_MATRICES

// One may define FCINT_SIMPLE_RECURSION to obtain the slow, basic recursion version for test purposes.
// Generally much slower, because no intermediate integrals of the recursion are ever stored in 
// lookup tables. The advantage of using this option (=defining the macro) are the very lean memeory requirements.
// In the language of the appendix in the article in J Phys Chem A cited above, not defining the macro implements
// a direct method, whereas defining it implements the "brute-force" method.
// #define FCINT_SIMPLE_RECURSION

#ifndef FCINT_USE_TRANSPOSED_MATRICES
#define FCINT_INDX(row,col,cols) ((col)+(row)*(cols))
#else
#define FCINT_INDX(row,col,cols) ((row)+(col)*(cols))
#endif
/**********************************************
 Inputs:
      N: the number of elements in the vector n
      Nprime: the number of elements in the vector n'
      n[]: the N dimensional vector with the quantum numbers <n1 n2 n3 n4...|
      nprime[]: the N' dimensional vector with the quantum numbers |n1' n2' n3' n4'...>
      RJdagg[][]: N by N' matrix
      twoR1[][]: N by N matrix
      RJD[]: N dimensional vector
      JR[][]: N' by N matrix
      twoP1[][]: N' by N' matrix
      onePD[]: N' dimensional vector
      OO: the overlap <00000|0'0'0'0'0'>
      init: nonzero if calling the FCint() with modified 'N,' 'Nprime', 'RJdagg', 'twoR1', 'RJD', 'JR', 'twoP1', 'onePD' or 'OO',
            zero if calling the FCint() with the same set of those variables (only some or all elements in 'n' or 'nprime' changed.
            The meaning of a nonzero 'init' is that this triggers removal of an internal database of previously calculated
            matrix elements. If 'init' is zero, previously calculated matrix elements are re-inserted into the recursion
            if possible (which is possible only if neither the dimensionality of the two vectors 'n' and 'nprime' nor
            any of the matrix or vector parameters was changed.)

 Return and output value: the overlap matrix <n1 n2 n3 n4....|n1' n2' n3' n4'>
************************************************/

#if ! defined FCINT_SIMPLE_RECURSION
// The definition of the integral, but without it's value or the associated matrices
// The values of N, Nprime are not needed within the class, but only convenient in braket_less.
class braket {
	private :
	public :
	vector<int> n ;
	vector<int> nprime ;
	braket(const int *nvect, const int dimn, const int *nprimevect, const int dimnp)
	{
		for(int i=0; i< dimn; i++)
			n.push_back(nvect[i]) ;
		for(int j=0; j< dimnp; j++)
			nprime.push_back(nprimevect[j]) ;
	}
} ;

// A way to sort instances of 'braket' in some way.
// We sort such that low values of n[] and nprime[] are kept early in the list, and whence found fast.
// Note that there is no internal check that the dimensions N and N' of the left or right are the same and
// that checks on bra1==bra2 and bra1!=bra2 are done by repeated calls of this operator with the two values swapped.
struct braket_less {
	bool operator()(const braket &bra1, const braket &bra2) const
	{
		int indx ;
		for(indx=0 ; indx < bra1.n.size() ; indx++) 
			if( bra1.n[indx] < bra2.n[indx])
				return true ;
			else if( bra1.n[indx] > bra2.n[indx])
				return false ;
		for(indx=0 ; indx < bra1.nprime.size() ; indx++) 
			if( bra1.nprime[indx] < bra2.nprime[indx])
				return true ;
			else if( bra1.nprime[indx] > bra2.nprime[indx])
				return false ;
		return false ;
	}
} ;

// A class that stores known (computed) integrals of type braket, including their values
class FCintheap {
	private :
	int N ;
	int Nprime ;
	map<braket,double,braket_less> oldFCs ;

	public :

	// default constructor
	FCintheap()
	{
		init() ;
	}

	// remove the old elements in the list of elements
	void init(int dimn=0, int dimnprime=0)
	{
		oldFCs.clear() ;
		N=dimn ;
		Nprime=dimnprime ;
	}

	// check for presence of a <bra|ket> in the list and return value in 'findVal' if found
	bool find(const int *n, const int *nprime, double *findVal)
	{
		braket tmp(n,N,nprime,Nprime) ;
		map<braket,double,braket_less>::iterator FCit = oldFCs.find(tmp) ;
		if( FCit != oldFCs.end() )	// found this braket in the list: store temporarily its value into 'findVal'
		{				// supposing that it's needed
			*findVal = FCit->second ;
			return true ;
		}
		else
			return false ;
	}

	// insert the <bra|ket> into the list of known values, using value 'insrtVal'
	void insert(const int *n, const int *nprime, const double insrtVal)
	{
		braket tmp(n,N,nprime,Nprime) ;
		oldFCs.insert(map<braket,double,braket_less>::value_type(tmp,insrtVal)) ;
	}
} ;


#endif

double FCint(const int *n, const int &N, const int *nprime, const int &Nprime, const double *RJdagg, const double *twoR1,
		const double *RJD, const double *JR, const double *twoP1, const double *onePD,const double &OO, const int &init)
{

#ifndef FCINT_SIMPLE_RECURSION
	static FCintheap knownFCs ;
#endif

#ifndef FCINT_SIMPLE_RECURSION
	// If requested, delete all older overlap integrals, and prepare to store NxNprime integrals
	if ( init)
		knownFCs.init(N,Nprime) ;
#endif

	// Search for <n,nprime> in the old list of integrals; if found, just return it
	double result ;
#ifndef FCINT_SIMPLE_RECURSION
	if( knownFCs.find(n,nprime,&result) )
		return result ;
#endif

	// At next deeper level of recursion, matrix elements can be re-used
	int nextinit=0 ;
	// search for non-zero n's in the bra with index i
	for( int i =0 ; i < N ; i++)
		if( n[i] )		// nonzero element in n found
		{
			int * newn= new int [N] ;
			int * newnk = new int [N] ;
			int * newnprime = new int [Nprime] ;
			// start with the negative last term in eq (8)
			// make a local copy of the n vector with the position 'i' decremented by 1, <n|
			memcpy(newn,n,N*sizeof(int)) ;
			newn[i]-- ;
				// recursive call to get <n|n'> in eq (8)
			result = -sqrt(2./n[i])*RJD[i]*
				FCint(newn,N,nprime,Nprime,RJdagg,twoR1,RJD,JR,twoP1,onePD,OO,nextinit) ;
				// first term in eq (8)
			for(int j=0; j< Nprime; j++)
	
				if( nprime[j] )		// test that n'-1_j is >= 0
				{
					memcpy(newnprime,nprime,Nprime*sizeof(int)) ;
					newnprime[j]-- ;		// subtract 1 at the j'th position, |n'-1_j
						// recursive call to get <n|n'-1_j> in eq (8)
					result += 2.*sqrt(nprime[j]/double(n[i]))
						*RJdagg[FCINT_INDX(i,j,Nprime)]
						*FCint(newn,N,newnprime,Nprime, RJdagg,twoR1,RJD,JR,twoP1,onePD,OO,nextinit) ;
				}
				// second term in eq (8)
			for(int k=0; k< N; k++)

				if( newn[k] )		// test that n-1_k stays >= 0
				{
					memcpy(newnk,newn,N*sizeof(int)) ;
					newnk[k]-- ;		// subtract 1 at the k'th position, <n-1_k|
						// recursive call to get <n-1_k|n'> in eq (8)
					result += sqrt(newn[k]/double(n[i]))
						*twoR1[FCINT_INDX(i,k,N)]
						*FCint(newnk,N,nprime,Nprime, RJdagg,twoR1,RJD,JR,twoP1,onePD,OO,nextinit) ;
				}
			
			delete [] newnprime ;
			delete [] newnk ;
			delete [] newn ;
#ifndef FCINT_SIMPLE_RECURSION
			knownFCs.insert(n,nprime,result) ;
#endif
			return result ;
		}
	// if we are here, all elements in n are already zero, so we search for a nonzero n' in the ket
	for( int j =0 ; j < Nprime ; j++)
		if( nprime[j] )		// nonzero element in n'
		{
			int * newn= new int [N] ;
			int * newnprime = new int [Nprime] ;
			int * newnprimek = new int [Nprime] ;
			// start with the last term in eq (9)
			// make a local copy of the n' vector with the position 'j' decremented by 1, |n'>
			memcpy(newnprime,nprime,Nprime*sizeof(int)) ;
			newnprime[j]-- ;
				// recursive call to get <n|n'> in eq (9)
			result = sqrt(2./nprime[j])*onePD[j]*
				FCint(n,N,newnprime,Nprime,RJdagg,twoR1,RJD,JR,twoP1,onePD,OO,nextinit) ;
				// first term in eq (9)
			for(int i=0; i< N; i++)
	
				if( n[i] )		// test that n-1_i is >= 0
				{
					memcpy(newn,n,N*sizeof(int)) ;
					newn[i]-- ;		// subtract 1 at the i'th position , <n-1_i|
						// recursive call to get <n-1_i|n'> in eq (9)
					result += 2.*sqrt(n[i]/double(nprime[j]))
						*JR[FCINT_INDX(j,i,N)]
						*FCint(newn,N,nprime,Nprime, RJdagg,twoR1,RJD,JR,twoP1,onePD,OO,nextinit) ;
				}
				// second term in eq (9)
			for(int k=0; k< Nprime; k++)

				if( newnprime[k] )		// test that n'-1_k stays >= 0
				{
					memcpy(newnprimek,newnprime,Nprime*sizeof(int)) ;
					newnprimek[k]-- ;		// subtract 1 at the k'th position
						// recursive call to get <n|n'-1_k> in eq (9)
					result += sqrt(newnprime[k]/double(nprime[j]))
						*twoP1[FCINT_INDX(j,k,Nprime)]
						*FCint(n,N,newnprimek,Nprime, RJdagg,twoR1,RJD,JR,twoP1,onePD,OO,nextinit) ;
				}
			
			delete [] newnprimek ;
			delete [] newn ;
			delete [] newnprime ;
#ifndef FCINT_SIMPLE_RECURSION
			knownFCs.insert(n,nprime,result) ;
#endif
			return result ;
		}
	// if we have come here, all n and nprime are zero, so we return the integral at the outermost leaves
	// of the recursion

#ifndef FCINT_SIMPLE_RECURSION
	knownFCs.insert(n,nprime,OO) ;
#endif
	return OO ;
}