/************************************************************************/
/*                                                                      */
/*  Program    :     B D K . c p p                                      */
/*                                                                      */
/*  Creation   :     30/05/02                                           */
/*                                                                      */
/*  Class      :     Main program                                       */
/*                                                                      */
/*  Programmer :     Uwe Trittmann                                      */
/*                                                                      */
/*  Purpose    :     PGR for calculating QCD_1+1 a la                   */
/*                   Bhanot/Demeterfi/Klebanov                          */
/************************************************************************/

#include<iostream>
using namespace std;
#include<cmath>
#include<string>
#include<vector>
#include<list>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<functional>
#include<time.h>
#include<limits.h>
#include<float.h>
#include "NAGdiagIV.h"          // diagonalization routines
//#include "diagonalization.h"

#define PI 3.141592653589793
#define ONE 2

#define KK 30              /* Maximal harmonic resolution */
#define max_nos 40000      /* maximal number of states */

int state[max_nos][KK+1];
int T_partner[max_nos];    // stores T partner of state and Z_2 QN
int nos;                   // number of states 
int nos_odd, nos_even;
int state_ket[max_nos][KK+3];
int matches[KK];
int indexHam,test;
int dim_even,dim_odd;      // dimensions of the Z_2 even/odd Hamiltonians
double *Ham; 	                       // variables for NAGLIB diagonalization 
double *eigenvalues,*eigenfunctions;   // Variables for diagonalization
double norm[max_nos];                  // norm of the states 

                
/*======================================================================== */

int mod(int i,int b)
{
  if (i>b) return i-b;
  else return i;
}

int perm(int n, int b)
{
  if (b%2!=0) return 1;
  else 
    {
      if (n%2!=0) return -1;
      else return 1;
    }     
}

int symmetry_factor(int n)
{
  int j,k,b,sym,equiv;

  sym = 0;  
  b = state[n][0];
  for (k=1; k<b; ++k) /* all possible permutations */
    {
      equiv = 1;
      j = 1;
      while ((equiv==1)&&(j<=b)) /* all bits until disagreement */
	{
	  if (state[n][mod(k+j,b)]!=state[n][j]) equiv = 0;
	  ++j;
	}
      if (equiv==1) ++sym;
    }
  return sym;
}

int skp(int out, int in)
/* Scalar product between state[out] and ket_state[in] */
/* returns 0 if SKP is 0 or # of total matches  */
{
  int j,k,b,sym,equiv;

  sym = 0;  
  b = state[out][0];
  for (k=0; k<b; ++k) /* all possible permutations */
    {
      equiv = 1;
      j = 1;
      while ((equiv==1)&&(j<=b)) /* all bits until disagreement */
	{
	  if (state_ket[in][mod(k+j,b)]!=state[out][j]) equiv = 0;
	  ++j;
	}
      if (equiv==1) 
	{
	  ++sym;
	  matches[sym-1]=k;
	}
    }
  return sym;
}

int test_existence(int nr)
/* tests if state vanishes due to fermionic statistics */
{
  if (state[nr][0]%2==1) return 1;
  if (fmod(1.0*state[nr][0]/(symmetry_factor(nr)+1),2.0)==0) return 1;
                                                                /* state OK */
  else return 0; /* state does not exist */
  
}

int test_state(int nr)
// test if state <nr> exists already in list of states in some permutation 
// return: 0 if not, (number of state + 1)*sign  if yes
{
  int i,j,k,found,equiv,b=state[nr][0];

  found = 0;
  i = 0;
  while ((found==0)&&(i<nos)) 
    {
      if (state[i][0]==b) 
	{
	  k=0;
	  while ((found==0)&&(k<b)) /* all possible permutations */
	    {
	      equiv = 1;
	      j = 1;
	      while ((equiv==1)&&(j<=b)) /* all bits until disagreement */
		{
		  if (state[i][mod(k+j,b)]!=state[nr][j]) equiv = 0;
		  ++j;
		}
	      if (equiv==1) found = 1;
	      ++k;
	    }
	}
      ++i;
    }
  if (found==0) return found;
  else 
    {
      if ((b%2==1)||((k-1)%2==0)) return i; // k cyclic rotations give sign 
      else return -i;                       // if b even and k odd
    }
}

void find_apBC_state(int b, int remaining_K)
// find states with odd momenta, i.e. consider fermions with anti-periodic 
// boundary conditions
{
  int n,m;

  for (n=1; n<=remaining_K; n+=2)
    {
      state[nos][b] = n;
      if (n<remaining_K) find_apBC_state(b+1,remaining_K-n);
      else 
	{
	  state[nos][0] = b;
	  if ((test_state(nos)==0) && (test_existence(nos)==1) && (b!=1))
	    {
	      ++nos;
	      for (m=1; m<=b-1; ++m) state[nos][m]=state[nos-1][m];
	      norm[nos-1]=1.0/sqrt(symmetry_factor(nos-1)+1.0);
	      if (test==1)
		{
		  printf("n=%3d: Tr [",nos-1);
		  for (m=1; m<=b; ++m) printf("B^+(%1d)",state[nos-1][m]);
		  printf("]/N^{%1d/2}",state[nos-1][0]);
		  if (symmetry_factor(nos-1)>0) 
		    printf("/sqrt{%1d}\n",symmetry_factor(nos-1)+1);
		  else printf("\n");
		}
	    }
	}
    }
}

void T_image(int nr)
// create image of state <nr> under T operation 
// (flip of color indices and sign), result is state <nos+1>
{
  int b=state[nr][0];

  state[nos+1][0]=b;
  state[nos+1][b]=state[nr][b];
  for (int i=1; i<b; ++i) state[nos+1][i]=state[nr][b-i];
}

int T_sign(int nr)
// find sign of state <nr> under T operation incl. fermion permutations
// return sign
{
  int b=state[nr][0],sign=1;

  for (int i=1; i<b-1; ++i)
    if ((b-i-1)%2==1) sign*=-1;
  if (b%2==1) sign*=-1;
  return sign;
}

void find_T_partners()
// find the partners of all states under T operation 
// (flip of color indices and sign)
{
  int sign,nr;

  dim_even=0;
  dim_odd=0;
  for (int i=0; i<nos; ++i)
    {
      T_image(i);
      nr=test_state(nos+1);
      sign=T_sign(i);
      if (nr<0) 
	{
	  sign*=-1;
	  nr*=-1;
	}
      //if (state[i][0]%2==1) sign*=-1;
      if ((nr-1)!=i)
	{
	  if ((nr-1)>i) 
	    {
	      T_partner[i]=nr;
	      ++dim_even;
	    }
	  else	    
	    {
	      T_partner[i]=-nr;
	      ++dim_odd;
	    }	
	}
      else
	{
	  if (sign>0)
	    {
	      T_partner[i]=nr;
	      ++dim_even;
	    }
	  else 
	    {
	      T_partner[i]=-nr;
	      ++dim_odd;
	    }
	}
      //printf("T(%2d)=%2d\n",i,T_partner[i]);
    }
  //printf("Even states: %2d\nOdd states:  %2d\n",dim_even,dim_odd);
}

void count_states()
{
  int i;

  nos_odd = nos_even = 0;
  for (i=0; i<=nos-1; ++i)
    {
      if (state[i][0]%2==0) ++nos_even;
      else ++nos_odd;
    }
  printf("nos_even= %3d\n",nos_even);
  printf("nos_odd = %3d\n",nos_odd);
}

double V_f(int n)
/* kinetic energy */
{
  int i,b;
  double sum;

  b = state[n][0];
  sum = 0.0;
  for (i=1; i<=b; ++i) sum += 1.0/state[n][i];
  return sum;
}

double T_f_1(int out, int n)
/* potential energy */
{
  int i,m,b;
  double sum;

  if (out!=n) return 0.0;
  b = state[n][0];
  sum = 0.0;
  for (i=1; i<=b; ++i) 
    for (m=1; m<=state[n][i]-2; m+=ONE) 
      sum += 2.0/(state[n][i]-m)/(state[n][i]-m);
  return sum;
}

double function_Tf2(int n1, int n2, int n4)
{
  if (n4-n2!=0) return 1.0/(n4-n2)/(n4-n2)-1.0/(n1+n2)/(n1+n2);
  else          return                    -1.0/(n1+n2)/(n1+n2);
}

double function_Tf3(int n1, int n2, int n3)
{
  return -1.0/(n3+n2)/(n3+n2)+1.0/(n1+n2)/(n1+n2);
}

double T_f_3(int out, int in)
/* potential energy, 3rd part */
/* B^+(k)B(l)B(i)B(j) */
{
  int i,j,b_in,b_out,nok,nosk,m,nr;
  int ket[100][KK+1];
  double sum,sum2;

  b_in  = state[in][0];
  b_out = state[out][0];
  sum = 0.0;
  if (b_in==b_out+2) sum = T_f_3(in,out); 
  if (b_in+2==b_out)
    {
      nok = 0;
      for (i=1; i<=b_in; ++i)
	{
	  ket[nok][b_in]=i;
	  ket[nok][0]=1;
	  if ((b_in%2==0)&&(i%2==0)) ket[nok][0] *= -1; 
	  for (j=1; j<b_in; ++j) ket[nok][j]=mod(i+j,b_in);
	  ++nok;
	}
      if (test==1)
	{
	  for (i=0; i<nok; ++i)
	    {
	      printf("|%2d> = |",i);
	      for (j=0; j<=b_in; ++j) printf("%2d ",ket[i][j]);
	      printf(">\n");
	    }
	}
      nosk = 0;
      for (nr=0; nr<nok; ++nr)
	for (i=1; i<=state[in][ket[nr][b_in]]-2; i+=ONE)
	  for (j=1; j<=state[in][ket[nr][b_in]]-i-1; j+=ONE)
	    {
	      state_ket[nosk][0]=ket[nr][0];
	      state_ket[nosk][1]=i;
	      state_ket[nosk][2]=j;
	      state_ket[nosk][3]=state[in][ket[nr][b_in]]-i-j;
	      state_ket[nosk][b_in+3]=state[in][ket[nr][b_in]]; 
	      for (m=4; m<=b_in+2; ++m)
		state_ket[nosk][m]=state[in][ket[nr][m-3]];
	      ++nosk;
	  }
      if (test==1)
	{
	  printf("state_out|%2d> = |",i);
	  for (j=0; j<=b_out; ++j) printf("%2d ",state[out][j]);
	  printf(">\n");
	  for (i=0; i<nosk; ++i)
	    {
	      printf("state|%2d> = |",i);
	      for (j=0; j<=b_in+3; ++j) printf("%2d ",state_ket[i][j]);
	      printf(">*%2d:",skp(out,i));
	      for (j=0; j<skp(out,i); ++j) printf("%2d",matches[j]);
	      printf("\n");
	    }
	}
      for (i=0; i<nosk; ++i) /* sum over all possible states */
	{
	  m = skp(out,i);
	  for (j=1; j<=m; ++j) /* sum over all matches */
	    {
	      sum2 = function_Tf3(state_ket[i][1],
				  state_ket[i][2],state_ket[i][3]);
	      if (test==1) printf("sum2 = %2d * %2d * %5.3f\n",
			       state_ket[i][0],perm(matches[j-1],b_in),sum2);
	      sum += state_ket[i][0]*perm(matches[j-1],b_in)*sum2;
	    }
	}
      sum *= norm[in]*norm[out];
    }
  if (test==1) printf("sum = %5.3f\n",sum);
  return sum;
}


double T_f_2(int out, int in)
/* potential energy, 2nd part, latest version */
/* B^+(k)^+B(l)B(i)B(j) */
{
  int i,j,nr,b_in,b_out,nok,nosk,m;
  int ket[100][KK];
  double sum,sum2;

  b_in  = state[in][0];
  b_out = state[out][0];
  sum = 0.0; 
  if (b_in==b_out)
    {
      nok = 0;
      for (i=1; i<=b_in; ++i)
	{
	  ket[nok][b_in-1]=i;
	  ket[nok][b_in]=mod(i+1,b_in);
	  ket[nok][0]=-1;
	  for (j=1; j<=b_in-2; ++j)
	    {
	      ket[nok][j]=mod(i+j+1,b_in);
	    }
	  if ((b_in%2==0)&&(i%2==0)) ket[nok][0]*=-1; 
	                      /* i= pos of 1+b_in-2, sign indicator, p.132 */
	  ++nok; 
	} 
      if (test==1)
	{
	  for (i=0; i<nok; ++i)
	    {
	      printf("|%2d> = |",i);
	      for (j=0; j<=b_in; ++j) printf("%2d ",ket[i][j]);
	      printf(">\n");
	    }
	}
      nosk = 0;
      for (nr=0; nr<nok; ++nr)
	for (i=1; i<=state[in][ket[nr][b_in-1]]+
	       state[in][ket[nr][b_in]]-1; i+=ONE)
	  {
	    state_ket[nosk][0]=ket[nr][0];
	    state_ket[nosk][1]=i;
	    state_ket[nosk][2]=state[in][ket[nr][b_in-1]]+
	      state[in][ket[nr][b_in]]-i;
	    state_ket[nosk][b_in+1]=state[in][ket[nr][b_in-1]]; 
	    /* add n1,n2 */
	    state_ket[nosk][b_in+2]=state[in][ket[nr][b_in]];
	    for (j=3; j<=b_in; ++j)
	      state_ket[nosk][j]=state[in][ket[nr][j-2]];
	    ++nosk;
	  }
      if (test==1)
	{
	  printf("state_out|%2d> = |",i);
	  for (j=0; j<=b_in; ++j) printf("%2d ",state[out][j]);
	  printf(">\n");
	  for (i=0; i<nosk; ++i)
	    {
	      printf("state|%2d> = |",i);
	      for (j=0; j<=b_in; ++j) printf("%2d ",state_ket[i][j]);
	      printf(">*%2d:",skp(out,i));
	      for (j=0; j<skp(out,i); ++j) printf("%2d",matches[j]);
	      printf("\n");
	    }
	}
      for (i=0; i<nosk; ++i) /* sum over all possible states */
	{
	  m = skp(out,i);
	  for (j=1; j<=m; ++j) /* sum over all matches */
	    {
	      sum2 = function_Tf2(state_ket[i][b_in+1],
				  state_ket[i][b_in+2],state_ket[i][2]);
	      if (test==1) printf("sum2 = %2d * %2d * %5.3f\n",
			      state_ket[i][0],perm(matches[j-1],b_in),sum2);
	      sum += state_ket[i][0]*perm(matches[j-1],b_in)*sum2;
	    }
	}
    }
  sum *= norm[in]*norm[out];
  if (test==1) printf("sum = %5.3f\n",sum);
  return sum;
}

void generate_Hamiltonian(int n, double x, int K, int Z, int test)
{
  int i,j,sign_in,sign_out;
  double norm_in,norm_out;

  indexHam=0;
  for (i=0; i<n; ++i) 
    {      
      if (Z*T_partner[i]>0)
	{
	  norm_in=1.0;
	  sign_in=1;
	  if (abs(T_partner[i])-1!=i)
	    { 
	      norm_in=1.0/sqrt(2.0);
	      T_image(i); // T image into <nos+1>
	      if (Z*test_state(nos+1)*T_sign(i)<0) sign_in=-1;  
	    }
	  for (j=0; j<n; ++j) 
	    if (Z*T_partner[j]>0)
	      {
		norm_out=1.0;
		sign_out=1;
		if (abs(T_partner[j])-1!=j) 
		  {
		    norm_out=1.0/sqrt(2.0);
		    T_image(j); // T image into <nos+1>
		    if (Z*test_state(nos+1)*T_sign(j)<0) sign_out=-1;  
		  }
		Ham[indexHam]=0.0;
		if (i==j) Ham[indexHam]+=V_f(i)*x;
		Ham[indexHam]+=2*(T_f_1(i,j)+T_f_2(i,j)+T_f_3(i,j));
		if (norm_in!=1.0)
		  {
		    if (abs(T_partner[i])-1==j) 
		      Ham[indexHam]+=sign_in*V_f(abs(T_partner[i])-1)*x;
		    Ham[indexHam]+=sign_in*2*(T_f_1(abs(T_partner[i])-1,j)
				      +T_f_2(abs(T_partner[i])-1,j)
				      +T_f_3(abs(T_partner[i])-1,j));
		    if (norm_out!=1.0)
		      {
			if (abs(T_partner[i])-1==abs(T_partner[j])-1) 
			  Ham[indexHam]+=
			       sign_in*sign_out*V_f(abs(T_partner[j])-1)*x;
			Ham[indexHam]+=sign_in*sign_out*2*(
			       T_f_1(abs(T_partner[i])-1,abs(T_partner[j])-1)
			      +T_f_2(abs(T_partner[i])-1,abs(T_partner[j])-1)
			      +T_f_3(abs(T_partner[i])-1,abs(T_partner[j])-1));
		      }
		  }
		if (norm_out!=1.0)
		  {
		    if (i==abs(T_partner[j])-1) 
		      Ham[indexHam]+=sign_out*V_f(abs(T_partner[j])-1)*x;
		    Ham[indexHam]+=sign_out*2*(T_f_1(i,abs(T_partner[j])-1)
				      +T_f_2(i,abs(T_partner[j])-1)
				      +T_f_3(i,abs(T_partner[j])-1));
		  }
		Ham[indexHam]*=norm_in*norm_out;
		++indexHam;
	      }
	}
    }
  if (test==1)
    {
      indexHam=0;
      for (i=0; i<n; ++i)
	{
	  if (Z*T_partner[i]>0)
	    {
	      for (j=0; j<n; ++j) 
		if (Z*T_partner[j]>0)
		  {
		    printf("%5.3f ",Ham[indexHam]);
		    ++indexHam;
		  }
	      printf("\n");
	    }
	}    
    }
}

void generate_unsym_Hamiltonian(int n, double x, int K, int test)
// generates Hamiltonian not Z_2 symmetrized
{
  int i,j;

  indexHam=0;
  for (i=0; i<n; ++i) 
    {
      for (j=0; j<n; ++j) 
	{  
	  Ham[indexHam]=0.0;
	  if (i==j) Ham[indexHam]=V_f(i)*x;
	  Ham[indexHam]+=2*(T_f_1(i,j)+T_f_2(i,j)+T_f_3(i,j));
	  ++indexHam;
	}
    }
  if (test==1)
    {
      indexHam=0;
      for (i=0; i<n; ++i)
	{
	  for (j=0; j<n; ++j) 
	    {
	      printf("%5.3f ",Ham[indexHam]);
	      ++indexHam;
	    }
	  printf("\n");
	}
    }
}

void display_states(int Z)
// display states in <Z> sector
{
  int m,b,count=0;

  for (int i=0; i<nos; ++i)
    if (Z*T_partner[i]>0)
      {
	b=state[i][0];
	printf("n=%3d: Tr[",count);
	for (m=1; m<=b; ++m) printf("B^+(%1d)",state[i][m]);
	printf("]/N^{%1d/2}",state[i][0]);
	if (norm[i]!=1.0) 
	  printf("/%5.3f\n",norm[i]+1);
	else printf("\n");
	if (abs(T_partner[i])-1!=i)
	  {
	    printf("         +Tr[");
	    for (m=1; m<=b; ++m) 
	      printf("B^+(%1d)",state[abs(T_partner[i])-1][m]);
	    printf("]/N^{%1d/2}",state[abs(T_partner[i])-1][0]);
	    if (norm[abs(T_partner[i])-1]!=1.0) 
	      printf("/%5.3f\n",norm[abs(T_partner[i])-1]+1);
	    else printf("\n");	
	  }
	++count;
      }
}

void default_test()
// tests if BDK 1993 results are reproduced
{
  int dim,i;

  printf("Default test ===============================================\n\n");
  printf("States and Hamiltonian (up to signs) should be the same \n");
  printf("as Eqs. (28)-(31) of BDK, Phys.Rev D48 (1993) 4980\n\n");
  test=0;
  printf("Z+ sector:\n-----------\n");
  find_apBC_state(1,10);
  find_T_partners();
  dim=dim_even;
  printf("%1d States:\n------------\n",dim);  
  display_states(1); 
  printf("\nZ+ Hamiltonian:\n---------------\n");
  test=0;
  Ham=dvector(0,dim*dim-1);
  generate_Hamiltonian(nos,0.0,10,1,1);
  printf("\n");
  diagonalize_symmetric(&Ham,dim,&eigenvalues,&eigenfunctions);
  printf("\n%1d Eigenvalues:\n---------------\n",dim);  
  for (i=0; i<dim; ++i) 
    {
      printf("EV[%2d]=%8.5f\n",i,10*eigenvalues[i]);
    }
  free_dvector(Ham,0,dim*dim-1);
  free_dvector(eigenvalues,0,dim-1);
  free_dvector(eigenfunctions,0,dim*dim-1);
  printf("\nZ- sector:\n----------\n");
  dim=dim_odd;
  printf("%1d States:\n------------\n",dim);  
  display_states(-1); 
  printf("\nZ- Hamiltonian:\n---------------\n");
  test=0;
  Ham=dvector(0,dim*dim-1);
  generate_Hamiltonian(nos,0.0,10,-1,1);
  printf("\n");
  diagonalize_symmetric(&Ham,dim,&eigenvalues,&eigenfunctions);
  printf("\n%1d Eigenvalues:\n---------------\n",dim);  
  for (i=0; i<dim; ++i) 
    {
      printf("EV[%2d]=%8.5f\n",i,10*eigenvalues[i]);
    }
  free_dvector(Ham,0,dim*dim-1);
  free_dvector(eigenvalues,0,dim-1);
  free_dvector(eigenfunctions,0,dim*dim-1);
  printf("\nEnd Default test ===========================================\n\n");
}

int main()
{
  int i,dim;
  int K,Z;
  double x;
  char name[100];
  FILE *fp;

  printf("\n\nQCD(1+1) with adjoint fermions [BDK, Phys.Rev D48 (1993) 4980]\n");
  printf("-------------------------------------------------------------------\n");
  printf("Anti-periodic boundary conditions!\n\n");

  default_test(); 
  
  K=0;
  printf("Harmonic resolution: K=");
  scanf("%d",&K);
  printf("Z_2 sector:          Z=");
  scanf("%d",&Z);
  printf("Mass parameter:      x=");
  scanf("%lf",&x);
  nos=0;
  test=0;
  find_apBC_state(1,K);
  find_T_partners();
  test=0;
  if (test==1) display_states(Z); 
  test=0;
  if (Z>0) dim=dim_even;
  else dim=dim_odd;
  printf("Dimension of Hamiltonian: %1d\n",dim);
  test=0;
  //dim=nos;
  Ham=dvector(0,dim*dim-1);
  generate_Hamiltonian(nos,1.0,K,Z,0);
  //generate_unsym_Hamiltonian(nos,1.0,K,0);
  diagonalize_symmetric(&Ham,dim,&eigenvalues,&eigenfunctions);
  if (Z>0) sprintf(name,"EVs_K%1dx%3.2fZp.dat",K,x);  
  else sprintf(name,"EVs_K%1dx%3.2fZm.dat",K,x);  
  fp=fopen(name,"w");
  printf("\nLowest eigenvalues (max. 100):\n------------------------------\n");
  for (i=0; i<dim; ++i) 
    {
      if (i<100) printf("EV[%1d]=%8.5f\n",i,K*eigenvalues[i]);
      fprintf(fp,"%1d %8.5f\n",K,K*eigenvalues[i]);
    }
  fclose(fp);
  printf("All eigenvalues in <%s>.\n",name);  
  if (Z>0) sprintf(name,"EFs_K%1dx%3.2fZp.dat",K,x);  
  else sprintf(name,"EFs_K%1dx%3.2fZm.dat",K,x);  
  printf("Store eigenfunctions? [neg.=>all, 0=>none] ");
  scanf("%d",&Z);
  if (Z!=0)
    {
      if ((Z<0)||(Z>dim)) Z=dim;
      fp=fopen(name,"w");
      fprintf(fp,"%1d ",Z);	  
      for (int j=0; j<Z; ++j) 
	fprintf(fp,"%8.5f ",K*eigenvalues[j]);
      fprintf(fp,"\n");	  
      for (i=0; i<dim; ++i) 
	{
	  fprintf(fp,"%1d",i);	  
	  for (int j=0; j<Z; ++j) 
	    fprintf(fp,"%8.5f ",eigenfunctions[i*dim+j]);
	  fprintf(fp,"\n");	  
	}      
      fclose(fp);      
      printf("Eigenfunctions in <%s>; eigenvalues in first row [Needs test!].\n",name);  
    }
  free_dvector(Ham,0,dim*dim-1);
  free_dvector(eigenvalues,0,dim-1);
  free_dvector(eigenfunctions,0,dim*dim-1);
  return 0;
}