#include <iostream>
#include <cstdlib>

#include <complex>
#include <vector>

#include "pSTO.h"

using namespace std ;

#define NMAX 10

#if 0
// copy constructor not needed because memberwise copy suffices
// Empty constructor without proper initialization. Do some minimum amount of work perhaps to recover from
// undetected hidden use by initialization by an invalid 'n'.
pSTO::pSTO() : Ylm()
{
	n= 999 ;
}
#endif

pSTO::pSTO(const int nn, const int ll, const int mm, const double expo, const complex<double> coef) 
		: Ylm(ll,mm,coef), n(nn), alpha(expo)
{
	if( nn >=0 && expo > 0. && abs(mm)<= ll)
	{
#ifdef USE_NNM
		if( n <= NMAX)
			Nnm = Nntab[n]*pow(2.*expo,n+0.5) ;
		else
		{
			cerr << "Unknown Nnm for pSTO n= " << nn << endl ;
			exit(EXIT_FAILURE) ;
		}
#endif
	}
	else
	{
		cerr << "Faulty initialization pSTO n= " << nn << " m= " << mm << endl ;
		exit(EXIT_FAILURE) ;
	}
}

pSTO & pSTO::operator*= (const double & fact)
{
	c *= fact ;
	return *this ;
}

pSTO & pSTO::operator*= (const complex<double> & fact)
{
	c *= fact ;
	return *this ;
}

#ifdef USE_NNM
/* sequence of 1/sqrt[(2n)!]
> interface(prettyprint=0) ;
> Digits := 30 ;
> for n from 0 to 11 do
>    evalf(1./( (2*n)!)^(1/2)) ;
> od ;
*/
const double pSTO::Nntab[NMAX+1] = { 
	1.,
	.707106781186547524400844362105,
	.204124145231931508183107006226,
	.372677996249964949401528944789e-1,
	.498011920555997349987007158205e-2,
	.524950656957260037797939633659e-3,
	.456910899277617353044394667089e-4,
	.338684891864543091458842524507e-5,
	.218620157633905874407878632395e-6,
	.124976825726157033253614529760e-7,
	.641117588536780424092550271881e-9
} ;
#endif
#undef NMAX

/** complex conjugation, including the coefficient
* Edmonds, (2.5.6)
* Y(l,-m) = (-1)^m *cc [Y(l,m)]
*/
pSTO conj(const pSTO & arg)
{
	pSTO resul(arg.n, arg.l, -arg.m, arg.alpha, std::conj(arg.c)) ;
	if ( abs(arg.m) %2 )	// (-1)^m=-1
		resul.c *= -1 ;
	return resul ;
}

// allow comparison operator to handle cases with slight numerical jitters in the exponents; otherwise the
// strict member-by-member comparison would work.  No need to compare Nnm or c.
bool operator== (const pSTO & first, const pSTO & secnd)
{
	if( first.n != secnd.n || Ylm(first) != Ylm(secnd) )
			return false ;
	// assuming DBL_DIG=15 digits (as in limits.h) of precision of a "double"
	if( fabs(first.alpha-secnd.alpha) > 1.e-14 * first.alpha)
		return false ;
	else
		return true ;
}
ostream & operator<<(ostream &os, const pSTO & someSTO)
{
	os << " [" << someSTO.n << " " << someSTO.l << " " << someSTO.m << ", " << someSTO.alpha << "] " ;
	return os ;
}

pSTO operator*(const complex<double> & fact, const pSTO & in)
{
	return pSTO(in.n, in.l, in.m, in.alpha, in.c*fact) ;
}

/** product of two STO's. The left one is *not* taken conj().
*/
vector<pSTO> operator*(const pSTO & left, const pSTO & right)
{
	const double gamma= left.alpha + right.alpha ;	// combined exponential
	const int n= left.n + right.n-1 ;		// combined new main r-power
	// cout << __LINE__ <<  " " << left << " " << right << endl;
	const vector<Ylm> ylm = Ylm(left)*Ylm(right) ;	// product expansion of the two spherical harmonics
	vector<pSTO> resul ;
	for(int i=0 ; i < (int) ylm.size() ; i++)
		resul.push_back( pSTO(n, ylm[i].l, ylm[i].m, gamma, ylm[i].c) ) ;
	return resul ;
}