#include<iostream>
#include<vector>
#include<string>
#include<math.h>
using namespace std;

struct data
{
	int n;
	vector<int>a;
	long double value;
	bool check;
};


// List of functions defined
long double bin(int n, int r);
long double bin2(int x, int y, int z);
long double N0(int k , vector<int>c); // This function will give number of rational curves in P_k ^2
int product(vector<int>a, vector<int>b); // This defines the intersection porduct of given two homology class in P_k ^2
vector< vector < int > > produceclock (int k, vector<int> c); 
vector< vector < int > > produceclock1 (int k, vector<int> c);
long double N1(int k , vector<int>beta); //This function will give number of elliptic curves in P_k ^2
long double M0(vector< int > beta); // This function will give number of rational curves in P1XP1
long double M1(vector< int > beta); //This function will give number of elliptic curves in P1XP1
vector<data>store10; // This function will store all the values of N0( beta)
vector<data>store11; // This function will store all the values of N1( beta)
vector<data> store00; // This function will store all the values of M0( beta)
vector<data> store01; // This function will store all the values of M1( beta)
vector<int> alpha;
/*
 int b;
 int c;
 int d;
*/
//***************************************************************************************

int main(){
    int k,n;
    cout.precision(30);
    string cont="y";
    while(cont == "Y" || cont =="y" || cont == "yes"){
//    b=0,c=0,d=0;
    cout<<"Enter the number of beta which you want to compute:"<<endl;
    cin>>n;
    while(n<= 0){
        cout<<"Kindly enter positive number."<<endl;
        cout<<"Enter the number of beta which you want to compute:"<<endl;
        cin>>n;
    }
    int check[n],a[n];
    vector< vector< int > > tempstorage;
    for(int p=0; p< n; p++)
    {  
        int choice=0;
        while(choice < 1 || choice > 2)
        {
            cout<<"enter 1 for P1XP1 and enter 2 for P^2 blown up upto 8 point for beta_"<<p+1<<endl;
            cin>> choice;
        }
        if(choice==1)
        {
            check[p]=1;
            cout<<"Kindly enter the homology class beta_"<<p+1<<" for P1XP1:"<<endl;
            vector<int> beta;
            int in;
            cin>>in;
            beta.push_back(in);
            cin>>in;
            beta.push_back(in);
            tempstorage.push_back(beta); 
        }
        if(choice==2)
        {
            check[p]=2;
            cout<<"Kindly enter the number of blowup points (call this number k). The maximum is 8."<<endl;
            cin>>a[p];
            cout<<"Specify the k+1 component integral vector beta_"<<p+1<<
            " (having only whole numbers). For example, 5 1 2 0 where 5 is the coefficient of L and the number of blowup points is 3."<<endl;
            int in;
            vector<int> beta;
            cin>>in;
            beta.push_back(in);
            for(int i=1; i< a[p]+1; i++) {
            cin>>in;
            beta.push_back(-in);
            }
            tempstorage.push_back(beta);
/*            
            	cout<<"c[i] is ";
			for(int i=0;i< beta.size() ;i++){
				cout<<beta[i]<<" ";
			}
			cout<<endl;
			vector < vector < int > > clock = produceclock1(beta.size()-1 ,beta);
			for(int i=0;i<clock.size();i++){
				cout<<clock[i][0]<<" ";
				for(int j=1;j<clock[i].size();j++){
					cout<<- clock[i][j]<<" ";
				}
				cout<<endl;
			}
*/			
        }
    }
    
    for(int p=0; p< n; p++)
    {
        if(check[p]==1)
        {
            cout<<"The number of rational curve of degree ("<< tempstorage[p][0]<<", "<<
            tempstorage[p][1]<<") in P1XP1 is"<<" "<<M0( tempstorage[p])<<endl;
            cout<<"The number of elliptic curve of degree ("<< tempstorage[p][0]<<", "<<
            tempstorage[p][1]<<") in P1XP1 is"<<" "<<M1( tempstorage[p])<<endl;
        }
        if(check[p]==2)
        {
            cout<<"The number of rational curve of degree ("<< tempstorage[p][0];
            if(a[p]!=0)
            {
                for(int i=1; i< a[p]+1; i++)
                cout<<", "<<-tempstorage[p][i];
            }
            cout<<") in P^2 blown up upto 8 points is ";
            cout<< N0(a[p],tempstorage[p])<< endl;
            alpha=tempstorage[p];
            cout<<"The number of elliptic curve of degree ("<< tempstorage[p][0];
            if(a[p]!=0)
            {
                for(int i=1; i< a[p]+1; i++)
                cout<<", "<<-tempstorage[p][i];
            }
            cout<<") in P^2 blown up upto 8 points is ";
            cout<<N1(a[p],tempstorage[p])<< endl;
            
        }
 /*       
        cout<<"Number of times Term 1 contribute is "<<b<<endl;
        cout<<"Number of times Term 2 contribute is "<<c<<endl;
        cout<<"Number of times Term 3 contribute is "<<d<<endl;
        while( ! alpha.empty())
        alpha.pop_back();
*/
    }
		    cout<<"Continue ? (y/n)"<<endl;
	        cin>>cont;
    }
        return 0;
}

// definitions of the functions
//****************************************************************
long double bin(int n, int r){
	
	if(r>n || r<0 || n<0) return 0;
	if(r==0 && n==0) return 1;
	 long double answer=1;
	for(int i=0;i<r;i++){
		answer=answer * (n-i);
	}
	long double fact=1;
	for(int i=0;i<r;i++){
		fact=fact*(i+1);
	}
	answer = answer/fact;
	 return answer;
}
//*********************************************************************

long double bin2(int x, int y, int z){
	
	return bin(x, y)*bin(x-y, z);
}
//*******************************************************************
vector< vector < int > > produceclock (int k, vector<int> c){
	vector< vector < int > > alphabetcprime;
	if(k==0) {
		for(int i=0;i<c[0]+1;i++){
				vector< int > temporary;
				temporary.push_back(i);
					alphabetcprime.push_back(temporary);
		}
		return alphabetcprime;
	}
	vector<int>smaller;
	for(int i=0;i<k;i++){
		smaller.push_back(c[i]);
	}
for(int i=0;i<c[k]+1;i++){
		vector< vector <int > > temporary = produceclock(k-1,smaller);
		//create vectors inductively
	for(int j=0;j<temporary.size();j++){
// Fix the last entry of the vector to be i
		temporary[j].push_back(i);
			}
	for(int j=0;j<temporary.size();j++){
		alphabetcprime.push_back(temporary[j]);
						}
}
return alphabetcprime;
}
//*********************************************************************
long double N0(int k , vector<int>c){
	for(int i=k; i > 0; i--)
	{
		if(c[i]==0)
		{
			c.pop_back();
			k--;
		}
		else break;
	}
	long double answer=0;
	long double tempsumi=0;
	data tempstorage;
	tempstorage.n=k;
	tempstorage.a= c;
	tempstorage.check= false;
		if(c[0]==0)
	{
	    for(int i=1;i<k+1;i++)
	    {
	        if(c[i]>0)
	        {
	            tempstorage.check=true;
	            tempstorage.value=0;
	            store10.push_back(tempstorage);
	            return tempstorage.value; 
	        }
	        tempsumi+=c[i];
	    } 
	        if(tempsumi==-1) 
	        {
	            tempstorage.check=true;
	            tempstorage.value=1;
	            store10.push_back(tempstorage);
	            return tempstorage.value; 
	        }
	        else
	        {
	            tempstorage.check=true;
	            tempstorage.value=0;
	            store10.push_back(tempstorage);
	            return tempstorage.value; 
	        }
	} // [0,...,1,0..] is 1 every other thing starting with 0 is 0.
	if(c[0]<0)
	{
	    tempstorage.check=true;
        tempstorage.value=0;
        store10.push_back(tempstorage);
        return tempstorage.value; 
	}  // [negative,...] is 0
	tempsumi=0;
	for(int i=1;i<k+1;i++) {
		if(c[i]<0) 
		{
		    tempstorage.check=true;
	        tempstorage.value=0;
	        store10.push_back(tempstorage);
	        return tempstorage.value; 
	        }
	} // [positive, ...., positive,...] is 0 So from now on c[0]>0 and the others leq 0
	if(c[0]==1) 
	{
	    for(int i=1;i<k+1;i++)
	    {
	        tempsumi+=c[i];
	    } 
	    if(tempsumi==0) 
	    {
	        tempstorage.check=true;
	        tempstorage.value=1;
	        store10.push_back(tempstorage);
	        return tempstorage.value; 
	        }
	} // [1,0,0,..] gives 1
	long double deltabeta = 3*c[0]-1;
	for(int u=1;u<k+1;u++){
									deltabeta -=c[u];									
								} 
if(deltabeta<0) 
{
    tempstorage.check=true;
    tempstorage.value=0;
    store10.push_back(tempstorage);
    return tempstorage.value; 
}
for(int r=0; r< store10.size();r++)
{
	if(store10[r].n== k && store10[r].a== c && store10[r].check== true)
	return store10[r].value;
}
if(deltabeta>=0 && deltabeta<3){
	vector< vector < int > > alphabetcprime;
	int loop=0;
	bool condi=true;
	while(condi==true){
				loop++;
		if(c[loop]>0){condi=false;}
	}
	vector<int>tci;
	for(int i=0;i<k+1;i++){
		if(i==loop) tci.push_back(c[loop]-1);
		if(i!=loop) tci.push_back(c[i]);
	}
alphabetcprime = produceclock(k, tci);
for(int i=0;i<alphabetcprime.size();i++){
	vector<int>tempone;
	vector<int>temptwo;
	tempone = alphabetcprime[i];
	for(int j=0;j<k+1;j++){
		temptwo.push_back(tci[j]-tempone[j]);
	}
	bool makesense=true;
							long double tempsum=0;
							long double anothertempsum=0;
							for(int j=0;j<k+1;j++){
								if(temptwo[j]<0 || tempone[j]>tci[j] || tempone[j]<0 || temptwo[j]>tci[j]) makesense=false;
								tempsum+=tempone[j];
								anothertempsum+=temptwo[j];
							}
							if(tempsum<=0 || anothertempsum<=0) makesense=false;
							long double deltabetaone = 3*tempone[0]-1;
							long double deltabetatwo = 3*temptwo[0]-1;
							for(int u=1;u<k+1;u++){
									deltabetaone -=tempone[u];
									deltabetatwo -=temptwo[u];
								} 
								if(deltabetaone<0 || deltabetatwo<0) makesense=false;
								if(tempone[0]<=0 || temptwo[0]<=0) makesense=false;
								if(makesense==true){
									long double betaonetwo = tempone[0]*temptwo[0];
								for(int u=1;u<k+1;u++){
									betaonetwo = betaonetwo - tempone[u]*temptwo[u];
								} 
								long double t=tempone[0]*(temptwo[0]*tempone[loop]*temptwo[loop]-tempone[0]*temptwo[loop]*temptwo[loop]);
								answer+=N0(k, tempone)*N0(k,temptwo)*betaonetwo*t*bin(deltabeta,deltabetaone);									
								}
	}
	answer+=(c[0]*c[0]-(c[loop]-1)*(c[loop]-1))*	N0(k,tci);
answer=answer/(c[0]*c[0] * c[loop]);
tempstorage.check=true;
tempstorage.value= answer;
store10.push_back(tempstorage);
return tempstorage.value;
}
if(deltabeta>=3){
vector< vector < int > > alphabetcprime;
alphabetcprime = produceclock(k, c);
for(int i=0;i<alphabetcprime.size();i++){
	vector<int>tempone;
	vector<int>temptwo;
	tempone = alphabetcprime[i];
	for(int j=0;j<k+1;j++){
		temptwo.push_back(c[j]-tempone[j]);
	}
	bool makesense=true;
							long double tempsum=0;
							long double anothertempsum=0;
							for(int j=0;j<k+1;j++){
								if(temptwo[j]<0 || tempone[j]>c[j] || tempone[j]<0 || temptwo[j]>c[j]) makesense=false;
								tempsum+=tempone[j];
								anothertempsum+=temptwo[j];
							}
							if(tempsum<=0 || anothertempsum<=0) makesense=false;
							long double deltabetaone = 3*tempone[0]-1;
							long double deltabetatwo = 3*temptwo[0]-1;
							for(int u=1;u<k+1;u++){
									deltabetaone -=tempone[u];
									deltabetatwo -=temptwo[u];
								} 
								if(deltabetaone<0 || deltabetatwo<0) makesense=false;
								if(tempone[0]<=0 || temptwo[0]<=0) makesense=false;
							if(makesense==true){		
								long double betaonetwo = tempone[0]*temptwo[0];
								long double Lbetaone = tempone[0];
								long double Lbetatwo = temptwo[0];
								for(int u=1;u<k+1;u++){
									betaonetwo = betaonetwo - tempone[u]*temptwo[u];
								}
								long double t = Lbetatwo*bin(deltabeta-3,deltabetaone-1)-Lbetaone*bin(deltabeta-3,deltabetaone);							
							 answer+=N0(k, tempone)*N0(k, temptwo)*betaonetwo*Lbetaone*t;
							}		
}
tempstorage.check=true;
tempstorage.value= answer;
store10.push_back(tempstorage);
return tempstorage.value;
}
}
//*********************************************************************

int product(vector<int>a, vector<int>b){
	long double answer= a[0]*b[0];
	if(a.size() <= b.size())
	{
		for(int i=1; i< a.size();i++)
		answer -= a[i]*b[i];
	}
	else
	{
		for(int i=1; i< b.size();i++)
		answer -= a[i]*b[i];
	}
	return answer;
}
//*********************************************************************
vector< vector < int > > produceclock1 (int k, vector<int> c){
	vector< vector < int > > alphabetcprime;
	vector<int> larger;
	larger.push_back(c[0]);
	for(int i=1; i<= k; i++)
	{
		larger.push_back(c[i]+1);
	}
	vector< vector <int> > temporary = produceclock(k, larger);
	for(int i= temporary.size()-1; i >=1; i--)
	alphabetcprime.push_back(temporary[i]);
	
	vector<int> positive;
	if(k!= 0)
	{
	positive.push_back(0);
	for(int i=1;i<=k; i++)
	positive.push_back(1);
	temporary= produceclock(k, positive);
	for(int i=1; i< temporary.size(); i++)
	{
		for(int j=1; j< temporary[i].size(); j++)
		{
			if(temporary[i][j] != 0)
		    temporary[i][j]=-1;
		}
		alphabetcprime.push_back(temporary[i]);
	}
}
	return alphabetcprime;
}
//*********************************************************************	
	
long double N1(int k , vector<int>beta){
	for(int i=k ; i > 0; i--)
	{
		if(beta[i] ==0)
		{
			beta.pop_back();
			k--;
		}
		else break;
	}
	vector<int> xi;
	xi.push_back(3);
	for(int i=1; i<=k; i++)
	{
		xi.push_back(1);
	}
	long double T1=0,T2=0,T3=0,T4=0, sum=0;
	int kappabeta = product(xi, beta);
	int dX= product(xi, xi);
	int b2X= 1+k;
//	cout<<dX<<endl;
	
	// initial conditions of the recursion
	
	if(beta[0] < 1 )
	{
		return 0;
	} 
	if(beta[0]> 1)
	{
		for(int i=1;i<=k;i++)
		{
			if(beta[i] < 0) // [positive, ...., positive,...] is 0 So from now on beta[0]>0 and the others leq 0
			{
		    return 0;
			}
			
		}
	}
	if(beta[0]==1)
	{
		long double temp=0;
		for(int i=1;i<=k;i++)
		temp += beta[i];
		
		if (temp==0 || temp == -1)
	        {
		    return 0;
			}
	}
	if(kappabeta < 2)
	{
		return 0;
	}
	
	data tempstorage;
	tempstorage.n=k;
	tempstorage.a= beta;
	tempstorage.check= false;
	
	for(int r=0; r< store11.size(); r++)
	{
		if(store11[r].n== k && store11[r].a== beta && store11[r].check== true)
		return store11[r].value;
	}
		vector< vector < int > > secondvector;
		vector< vector < int > > thirdvector;
		secondvector= produceclock1(k, beta);
		thirdvector= produceclock1(k, beta);
		
		for(int i=0; i< secondvector.size(); i++)
		{
			// computation of T1
		for(int j=0; j< thirdvector.size(); j++)
		{
			vector<int> beta11;
			vector<int> beta12;
			vector<int> beta13;
			bool makesense1 = true;
			long double t1=0;
			int count1=0;
			beta12= secondvector[i];
			beta13= thirdvector[j];
			for(int p=0; p<k+1; p++)
			beta11.push_back(beta[p]- beta12[p]- beta13[p] );
			int kappabeta11 = product(xi, beta11);
			int kappabeta12 = product(xi, beta12);
			int kappabeta13 = product(xi, beta13);
			
			// Put some restriction on beta1,beta2,beta3
			if(beta11[0] <= 0 )
				{
					makesense1 = false;
				}
			// Next condition is beta13 not equal to [0,0,0.....,0]
			for(int i=0; i< beta11.size();i++)
				count1 += beta11[i]* beta11[i];
			if(count1==0) 
			makesense1 = false;
			if(makesense1== true)
			{
				t1 = 12* bin2(kappabeta-2, kappabeta12-1, kappabeta13-1)* 2* kappabeta12* pow(kappabeta13,2)* product(beta11,beta12)*
				(( ((4* kappabeta11) + kappabeta12 - (2* kappabeta13) ) * product(beta12,beta13) )- (3* kappabeta12 * product(beta11,beta13)))* N0(k,beta12)* N0(k,beta13);
				if(t1!= 0)
				t1= t1* N1(k, beta11);
				T1 += t1;
/*				
				bool print1=true;
				for(int i=0; i< k+1; i++)
				{
				    if(beta11[i]+ beta12[i]+beta13[i] != alpha[i])
				    {
				       print1=false;
				       break;
				    }
				}
				if(t1 != 0 && print1 == true)
				{
				    b++;
					cout<<"Term 1 is \n";
					cout<<"beta1 is (";
					cout<<beta11[0];
					if(k != 0)
					{
						cout<<", ";
						for(int i=1; i<k;i++)
				    	cout<<-beta11[i]<<", ";
				    	cout<<-beta11[k];
					}
					cout<<")"<<endl;
					cout<<"beta2 is (";
					cout<<beta12[0];
					if(k != 0)
					{
					cout<<", ";	
					for(int i=1; i<k;i++)
					cout<<-beta12[i]<<", ";
					cout<<-beta12[k];
					}
					cout<<")"<<endl;
					cout<<"beta3 is (";
					cout<<beta13[0];
					if(k != 0)
					{
						cout<<", ";
					    for(int i=1; i<k;i++)
					    cout<<-beta13[i]<<", ";
				    	cout<<-beta13[k];
					}
					cout<<")"<<endl;
					cout<<"value= "<< t1/12<<" %%%%%%%%%%% "<<endl;
				}
*/				
			}
		}
		   // computation of T2
		    vector<int> beta21;
			vector<int> beta22;
			long double t2=0;
			int count2=0;
			bool makesense2= true;
			beta22= secondvector[i];
			for(int j=0; j<k+1; j++)
			beta21.push_back(beta[j]- beta22[j] );
			int kappabeta21 = product(xi, beta21);
			int kappabeta22 = product(xi, beta22);
			// Put some restriction on beta1,beta2
			if(beta21[0] < 0 ) //beta21[0] can't be negative
				{
					makesense2 = false;
				}
			// Next condition is beta22 not equal to [0,0,0.....,0]
			for(int l=0; l< beta22.size();l++)
				count2 += beta22[l]* beta22[l];
				if(count2 == 0) 
				{
					makesense2 = false;
				}
				count2=0;
					
			if(makesense2== true)
			{
				t2= 12* ((bin(kappabeta -2, kappabeta21-1) *4* pow(kappabeta22,2)*((2*kappabeta21* kappabeta22)- pow(kappabeta22,2)- (3* dX * product(beta21,beta22))))+
				(bin(kappabeta-2,kappabeta21)* ((2* product(beta21,beta22)* kappabeta22 * dX * ((4* kappabeta21)+ kappabeta22)) + (4 * kappabeta21* pow(kappabeta22,2)*
				((2* kappabeta21)- kappabeta22)) ))) * N0(k, beta22); 
				if(t2 !=0)
				{
					t2= t2* N1(k, beta21); // this is the formula for T2
					if(t2 != 0)
					{
						T2 += t2;
					}
				}
/*				
				bool print2=true;
				for(int i=0; i<k+1; i++)
				{
				    if(beta21[i] + beta22[i] != alpha[i])
				    {
				        print2= false;
				        break;
				    }
				}
				if(t2 != 0 && print2 ==true)
				{
				    c++;
					cout<<"Term 2 is \n";
					cout<<"beta1 is (";
					cout<<beta21[0];
					if(k!= 0)
					{
						cout<<", ";
						for(int i=1; i<k;i++)
					    cout<<-beta21[i]<<", ";
				    	cout<<-beta21[k];
					}
					cout<<")"<<endl;
					cout<<"beta2 is (";
					cout<<beta22[0];
					if(k!=0)
					{
						cout<<", ";
						for(int i=1; i<k;i++)
					    cout<<-beta22[i]<<", ";
					    cout<<-beta22[k];
					}
				
					cout<<")"<<endl;
					cout<<"value= "<<t2/12<<" $$$$$$$$$$$ "<<endl;
				}
*/
			}
			//computation of T3
			vector<int> beta31;
			vector<int> beta32;
			long double t3=0;
			int count3=0;
			bool makesense3= true;
			beta32= secondvector[i];
			for(int p=0; p<k+1; p++)
			beta31.push_back(beta[p]- beta32[p]);
/*
			for(int p=0; p< k+1;p++)
			{
			    beta31.push_back(beta21[p]);
			    beta32.push_back(beta22[p]);
			}
*/			
			int kappabeta31= product(xi, beta31);
			int kappabeta32= product(xi, beta32);
			// Put some restriction on beta1,beta2
			if(beta31[0] < 0)
				{
					makesense3 = false;
				}
			// Next condition is beta32 not equal to [0,0,0.....,0]	
			for(int t=0; t< beta32.size();t++)
			count3 += beta32[t]* beta32[t];
			if(count3 == 0) 
				makesense3 = false;
				count3=0;
			if(makesense3== true)
			{
				t3= bin(kappabeta-2, kappabeta31-1)* pow(kappabeta32,2)* product(beta31,beta32)*(pow(kappabeta31,2)*(kappabeta31- (2* kappabeta32)- (6*product(beta31, beta32)))+
				(kappabeta32* product(beta31,beta31)*((4*kappabeta31)+ kappabeta32)))* N0(k,beta31) * N0(k,beta32); //This is the formula for T3		
				if(t3 != 0)
				  {
//					cout<<"t3 before counting "<<-t3<<endl;
//					cout<<"partial term before counting "<< -T3<<endl;
					T3 = T3+ t3;
//					cout<<"t3 after counting "<<-t3<<endl;
//					cout<<"partial term after counting "<< -T3<<endl;
                   }
/*                   
				bool print3=true;
				if(beta.size() == alpha.size() )
				{
				for(int t=0; t< k+1;t++)
				{
					if(beta31[t]+beta32[t] != alpha[t])
					{
						print3= false;
						break;
					}
				}
			    }
			    if(t3 != 0 && print3==true)
				{					
					d++;
					cout<<"Term 3 is \n";
					cout<<"beta1 is (";
					cout<<beta31[0];
					if(k!=0)
					{
						cout<<", ";
						for(int t=1; t<k;t++)
			    		cout<<-beta31[t]<<", ";
			    		cout<<-beta31[k];
					}
					cout<<")"<<endl;
					cout<<"beta2 is (";
					cout<<beta32[0];
					if(k!=0)
					{
						cout<<", ";
						for(int t=1; t<k;t++)
				    	cout<<-beta32[t]<<", ";
				    	cout<<-beta32[k];
					}
					cout<<")"<<endl;					
					cout<<" value= "<<-t3<<" ********* "<<endl;
				}
*/						
			}
		}
		//computation of T4
		T4= pow(kappabeta, 3)* (((2 + b2X)* kappabeta) - dX) * N0(k, beta);
/*		
		bool print34=true;
		if(beta.size() == alpha.size())
		{
		for(int i=0; i<k+1; i++)
		{
		    if(beta[i] != alpha[i])
		    print34=false;
		    break;
		}
	}
		if(T3+T4 !=0 && print34 == true){
		cout<<"Sum of Term 3 and Term 4 is "<<- (T3+T4)/12<<endl;
	}
	if(print34==true && beta.size() == beta.size() )
	{
		cout<<"Total contribution of term 1 is "<<T1/12<<endl;
		cout<<"Total contribution of term2 is "<<T2/12<<endl;
		cout<<"Total contribution of term3 is "<<-T3<<endl;
		cout<<"Total contribution of term4 is "<<-T4<<endl;
     }
*/     
		sum= (T1 + T2 - T3 -T4)/ (6*12* pow(dX, 2));
		tempstorage.value =sum;
		tempstorage.check= true;
		store11.push_back(tempstorage);
		return tempstorage.value;
}
//*********************************************************************
long double M0(vector< int > beta)
{
	long double sum=0;
	long double kappabeta= 2*(beta[0] + beta[1]);
	data tempstorage;
	tempstorage.n=1;
	tempstorage.check=false;
	// initial codition of the recurtion
	if((beta[0]==1 && beta[1] == 0) || (beta[0]==0 && beta[1] == 1))
	{
	    tempstorage.check=true;
	    tempstorage.value=1;
	    tempstorage.a= beta;
	    store00.push_back(tempstorage);
	    return tempstorage.value;
	}
	if(beta[0] < beta[1])
	{
	    vector<int> betaprime;
	    betaprime.push_back(beta[1]);
        betaprime.push_back(beta[0]);
	    tempstorage.a= betaprime;
	    tempstorage.check= true;
	    tempstorage.value = M0(betaprime);
	    store00.push_back(tempstorage);
	}
	for(int k=0; k< store00.size(); k++)
	{
	    if(store00[k].a== beta && store00[k].check==true && store00[k].n==1)
	    return store00[k].value;
	}
	vector< vector <int> > firstvector= produceclock(1, beta);
	for(int i=1; i< firstvector.size()-1; i++)
	{
	    vector< int > beta1;
    	vector< int > beta2;
    	beta1= firstvector[i];
    	beta2.push_back(beta[0]- beta1[0]);
    	beta2.push_back(beta[1]- beta1[1]);
    	long double kappabeta1= 2*(beta1[0] + beta1[1]);
    	long double kappabeta2= 2*(beta2[0] + beta2[1]);
    	long double beta1dotbeta2 = (beta1[0]*beta2[1]) + (beta2[0]* beta1[1]);
    	long double t=0;
    	t= beta1dotbeta2* beta2[1]*((beta1[0]* bin(kappabeta-4,kappabeta1-2))- (beta2[0]* bin(kappabeta-4, kappabeta2-1)))*
    	M0(beta1)* M0(beta2);
    	sum += t;
	}
	tempstorage.check=true;
	tempstorage.value= sum;
	tempstorage.a= beta;
	store00.push_back(tempstorage);
	return tempstorage.value;
}
//*********************************************************************

long double M1(vector< int > beta)
{
    long double sum=0,dX=8,b2X=2, T1=0,T2=0,T3=0,T4=0;
	long double kappabeta= 2*(beta[0] + beta[1]);
	data tempstorage;
	tempstorage.n=1;
	tempstorage.check=false;
	// initial codition of the recurtion
	if((beta[0]==1 && beta[1] == 0) || (beta[0]==0 && beta[1] == 1))
	{
	    tempstorage.check=true;
	    tempstorage.value=0;
	    tempstorage.a= beta;
	    store01.push_back(tempstorage);
	    return tempstorage.value;
	}
	if(beta[0] < beta[1])
	{
	    vector<int> betaprime;
	    betaprime.push_back(beta[1]);
        betaprime.push_back(beta[0]);
	    tempstorage.a= betaprime;
	    tempstorage.check= true;
	    tempstorage.value = M1(betaprime);
	    store01.push_back(tempstorage);
	}
	for(int k=0; k< store01.size(); k++)
	{
	    if(store01[k].a== beta && store01[k].check==true && store01[k].n==1)
	    return store01[k].value;
	}
		vector< vector < int > > secondvector;
		vector< vector < int > > thirdvector;
		secondvector= produceclock(1, beta);
		thirdvector= produceclock(1, beta);
		
		for(int i=1; i< secondvector.size()-1; i++)
		{
		    //computation of T1
		    for(int j=1; j< thirdvector.size()-1; j++)
		    {
		        vector<int> beta11;
			    vector<int> beta12;
			    vector<int> beta13;
		    	bool makesense1 = true;
		    	long double t1=0;
		    	beta12= secondvector[i];
		    	beta13= thirdvector[j];
		    	beta11.push_back(beta[0]- beta12[0]- beta13[0]);
		    	beta11.push_back(beta[1]- beta12[1]- beta13[1]);
			// Put some restriction on beta1,beta2,beta3
			if((beta11[0] < 0 || beta11[1] < 0) || (beta11[0]==0 && beta11[1]==0))
				{
					makesense1 = false;
				}
			if(makesense1== true)
			{
			    long double kappabeta11 = 2*(beta11[0] + beta11[1]);
		    	long double kappabeta12 = 2*(beta12[0] + beta12[1]);
		    	long double kappabeta13 = 2*(beta13[0] + beta13[1]);
		    	long double beta11dotbeta12 = (beta11[0]*beta12[1]) + (beta12[0]* beta11[1]);
		    	long double beta11dotbeta13 = (beta11[0]*beta13[1]) + (beta13[0]* beta11[1]);
		    	long double beta12dotbeta13 = (beta12[0]*beta13[1]) + (beta13[0]* beta12[1]);
		    	
				t1 = 12* bin2(kappabeta-2, kappabeta12-1, kappabeta13-1)* 2* kappabeta12* pow(kappabeta13,2)* beta11dotbeta12*
				(( ((4* kappabeta11) + kappabeta12 - (2* kappabeta13) ) * beta12dotbeta13 )- (3* kappabeta12 * beta11dotbeta13))* M0(beta12)* M0(beta13);
				if(t1!= 0)
				t1= t1* M1(beta11);
				T1 += t1;
/*				
				if(t1 != 0)
				{
					cout<<"Term 1 is \n";
                    cout<<"beta1 is ("<<beta11[0]<<", "<<beta11[1]<<")"<<endl;
                    cout<<"beta2 is ("<<beta12[0]<<", "<<beta12[1]<<")"<<endl;
                    cout<<"beta1 is ("<<beta13[0]<<", "<<beta13[1]<<")"<<endl;
					cout<<" value= "<<t1/12<<" %%%%%%%%%%% "<<endl;
				}
*/
			}
		    }
		    // computation of T2
		    vector<int> beta21;
			vector<int> beta22;
			long double t2=0;
			bool makesense2= true;
			beta22= secondvector[i];
			beta21.push_back(beta[0]- beta22[0] );
			beta21.push_back(beta[1]- beta22[1] );
			long double kappabeta21 = 2*(beta21[0] + beta21[1]);
			long double kappabeta22 = 2*(beta22[0] + beta22[1]);
			// Put some restriction on beta1,beta2
			if(beta21[0] < 0 || beta21[1] < 0 || (beta21[0]==0 && beta21[1]==0))
				{
					makesense2 = false;
				}
				if(makesense2==true)
				{
			        long double beta21dotbeta22 = (beta21[0]*beta22[1]) + (beta22[0]* beta21[1]);
			        
				    t2= 12* ((bin(kappabeta -2, kappabeta21-1) *4* pow(kappabeta22,2)*((2*kappabeta21* kappabeta22)- pow(kappabeta22,2)- (3* dX * beta21dotbeta22)))+
				    (bin(kappabeta-2,kappabeta21)* ((2* beta21dotbeta22 * kappabeta22 * dX * ((4* kappabeta21)+ kappabeta22)) + (4 * kappabeta21* pow(kappabeta22,2)*
				    ((2* kappabeta21)- kappabeta22)) ))) * M0(beta22); 
				   if(t2 !=0)
			       	{
					t2= t2* M1(beta21); // this is the formula for T2
					if(t2 != 0)
					{
						T2 += t2;
					}
				}
/*				
				if(t2 != 0)
				{
					cout<<"Term 2 is \n";
					cout<<"beta1 is ("<<beta21[0]<<", "<<beta21[1]<<")"<<endl;
					cout<<"beta2 is ("<<beta22[0]<<", "<<beta22[1]<<")"<<endl;
					cout<<" value= "<<t2/12<<" $$$$$$$$$$$ "<<endl;
				}
*/
			}
			//computation of T3
			vector<int> beta31 = beta21;
			vector<int> beta32 = beta22;
			long double t3=0;
			bool makesense3= true;
/*			
            beta32= secondvector[i];
			for(int p=0;p<2;p++)
			beta31.push_back(beta[p]- beta32[p]);
*/
			// Put some restriction on beta1,beta2
			if((beta31[0] < 0 || beta31[1] < 0) || (beta31[0]==0 && beta31[1]==0))
				{
					makesense3 = false;
				}
			if(makesense3== true)
			{
			    long double kappabeta31= kappabeta21;
			    long double kappabeta32= kappabeta22;
			    long double beta31dotbeta32 = (beta31[0]*beta32[1]) + (beta32[0]* beta31[1]);
			    long double beta31dotbeta31 = 2* beta31[0]*beta31[1];
				t3= bin(kappabeta-2, kappabeta31-1)* pow(kappabeta32,2)* beta31dotbeta32*(pow(kappabeta31,2)*(kappabeta31- (2* kappabeta32)- (6*beta31dotbeta32))+
				(kappabeta32* beta31dotbeta31*((4*kappabeta31)+ kappabeta32)))* M0(beta31) * M0(beta32); //This is the formula for T3
				T3 += t3;
/*				
			    if(t3 != 0)
				{
					cout<<"Term 3 is \n";
					cout<<"beta1 is ("<<beta31[0]<<", "<<beta31[1]<<")"<<endl;
					cout<<"beta2 is ("<<beta32[0]<<", "<<beta32[1]<<")"<<endl;
					cout<<" value= "<<-t3/12<<" ********* "<<endl;
				}
*/

			}
		}
		//computation of T4
		T4= pow(kappabeta, 3)* (((2 + b2X)* kappabeta) - dX) * M0(beta);
/*		
		if(T4!= 0)
		{
			cout<<"Term 4 is "<<- T4/12<<endl;
		}
*/		
		sum= (T1 + T2 - T3 -T4)/ (6*12* pow(dX,2));
		tempstorage.value =sum;
		tempstorage.check= true;
		tempstorage.a= beta;
		store01.push_back(tempstorage);
		return tempstorage.value;
}