//////////////////////////////////////////////////////////////////////////
//                                                                      //
//  Program    :     S Y M _ w i t h _ a p B C . c p p                  //
//                                                                      //
//  Creation   :     April 2000, revised August 2001                    //
//                                                                      //
//  Class      :     Main program                                       //
//                                                                      //
//  Programmer :     Uwe Trittmann                                      //
//                                                                      //
//  Purpose    :     PGR for calculating the spectrum of SYM(1+1) with  //
//                   anti-periodic boundary conditions, see paper       //
//                   S.Pinsky, U.Trittmann, Phys.Rev. D62 (2000) 087701 //
//                   (hep-th/0005055).                                  //
//                                                                      //
//----------------------------------------------------------------------//
//                                                                      //
// Subroutines :     NAGdiagII.h      ; diagonalization routines        //
//                   allocation.h     ; allocation of vectors           //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include <malloc.h>
#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include <iostream>
#include "NAGdiagII.h"          // diagonalization routines
 
#define KKK 10                  // maximal KK 
#define max_nos 15000
#define ef_relevance 1e-3       // cutoff for showing ef contribution 
#define TRUNC 10                // parton number trunc., TRUNC=KKK: no trunc.
#define TruncType 0             // 0/1 bosons/fermions don't make sense 

int KK;                         // harmonic resolution (must be even!) 
int state[max_nos][TRUNC+4];    // the states 
int out_state[max_nos][TRUNC+4];// out-states 
double amplitude[max_nos];      // amplitude of out-states 
double norm[max_nos];           // norm of states 
int nos,nos2,nos3,nos4;         // number of out-states 
int sector[TRUNC+4];            // array of sector boundarys  
int K_min;                      // minimal K, K_min >= 4 
int noos;                       // number of out-states 
int off[TRUNC+4];               // unused momenta for permutations 
int R[TRUNC+4];                 // return of function R() 
int perm_nr;                    // permutation number for get_state() 
double *Ham,*Ham2,*Ham3;        // variables for NAGLIB diagonalization 
double *r,*v,*e;                // variables for NAGLIB diagonalization 
int test;                       // test==1 gives print-outs 
int *numb;                      // ef number after sorting of evs 

//--------------------------------------------------------------------

void free_Intvector(int* v, long nl, long nh)
// free a int vector allocated with intvector() 
{
	free((char*) (v+nl-1));
}

void init_Ham(int dim)
// initialize Hamiltonian matrix for NAGLib diagonalization 
{  
	Ham = dvector(0,dim*dim - 1);
	Ham2= dvector(0,dim*dim - 1);
	Ham3= dvector(0,dim*dim - 1);
}

void free_Ham(int dim)
// initialize Hamiltonian matrix for NAGLib diagonalization 
{  
  free_dvector(Ham,0,dim*dim - 1);
  free_dvector(Ham2,0,dim*dim - 1);
  free_dvector(Ham3,0,dim*dim - 1);
}

// ----------------------------------------------------------------------

int mod(int i, int b)
// calculate i modulo b 
{
  if (i<0) mod(i+b,b);
  if (i>b) return mod(i-b,b);
  else return i;
}


int RRR(int j, int n, int nr)
// calculate shifted momentum string, input in <off[]>, output in <R[]> 
// <n>-<off[0]> momenta, <j> shifts to right 
{
  int i,k,count=0,found;
  
  for (i=1; i<=n; ++i)
    {      
      found=0;
      for (k=1; k<=off[0]; ++k) 
	if (off[k]==i+j) 
	  {
	    ++found; 
	    break;
	  }
      if (found==0)
	{
	  /*R[mod(count-j+state[nr][0]-off[0],state[nr][0]-off[0])-1]
	    =state[nr][i];*/
	  R[count]=state[nr][mod(i+j,state[nr][0])];
	  ++count;
	}
    }
}

int symmetry_factor(int n)
/* calculate symmetry factor of state, no symm. => s=0 */
{
  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 test_existence(int nr)
/* tests if state vanishes due to fermionic statistics */
{
  if ((state[nr][0]%2!=0)||(state[nr][0]/(symmetry_factor(nr)+1))%2==0) 
    return 1;    /* state OK */
  else return 0; /* state does not exist */
}

int test_state(int b)
{
  int i,j,k,found,equiv;

  found = 0;
  i = 0;
  while ((found==0)&&(i<nos)) 
    {
      if (state[i][0]==b) 
	{
	  k=1;
	  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[nos][j]) equiv = 0;
		  ++j;
		}
	      if (equiv==1) found = 1;
	      ++k;
	    }
	}
      ++i;
    }
  return found;
}


void find_state(int b, int remaining_K)
/* find states */
{
  int n,m;
  double sym;

  for (n=1; n<=remaining_K; n+=2)
    {
      state[nos][b] = n;
      if (n<remaining_K) find_state(b+1,remaining_K-n);
      else 
	{
	  state[nos][0] = b;
	  if ((test_state(b)==0) && (b!=1) &&  test_existence(nos))
	    {
	      ++nos;
	      for (m=1; m<=b-1; ++m) state[nos][m]=state[nos-1][m];	      
	      if (nos%100==0) printf("|%1d>\n",nos);

	      sym=symmetry_factor(nos-1);
	      if (sym>0) norm[nos-1]=1.0/sqrt(sym+1.0);
	      else norm[nos-1]=1.0;	      
	    }
	}
    }
}


void find_state_trunc(int b, int remaining_K, int p)
/* find states with up to <p> partons */
{
  int n,m;
  double sym;

  for (n=1; n<=remaining_K; n+=2)
    {
      state[nos][b] = n;
      if ((n<remaining_K)&&(b<=p)) find_state_trunc(b+1,remaining_K-n,p);
      else 
	{
	  state[nos][0] = b;
	  if ((test_state(b)==0) && (b!=1) && test_existence(nos) && (b<=p))
	    {
	      ++nos;
	      for (m=1; m<=b-1; ++m) state[nos][m]=state[nos-1][m];
	
	      if (nos%100==0) printf("|%1d>\n",nos);

	      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]);
		}
	      sym=symmetry_factor(nos-1);
	      if (sym>0) norm[nos-1]=1.0/sqrt(sym+1.0);
	      else norm[nos-1]=1.0;
	      
	      if (test==1)
		{
		  if (symmetry_factor(nos-1)>0) 
		    printf("/%5.3f\n",norm[nos-1]);
		  else printf("\n");
		}
	    }
	}
    }
}


int get_state(int nr)
/* find state that is cyclic permutation of out_state[nr] */
{
  int i,j,k,b,found,equiv;

  found=0;
  i=0;
  b=out_state[nr][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)]!=out_state[nr][j]) equiv = 0;
		  ++j;
		}
	      if (equiv==1) found = 1;
	      ++k;
	    }
	}
      ++i;
    }
  if (!found) return -1;
  else 
    {
      perm_nr=b-k+1;
      if (perm_nr==b) perm_nr=0;
      return i-1;
    }
}


void Q_minus(int plus, int nr)
/* act with Q^-<plus>1/2 on state <nr>, <plus>=+/-1 */
{
  int i,k,j;
  
  for (i=1; i<=state[nr][0]; ++i)
    {
      for (k=1; k<state[nr][i]+plus; k+=2)
	{
	  amplitude[noos]=1;
	  if ((i+1)%2!=0) amplitude[noos]*=-1;
	  if ((state[nr][0]*(i+1))%2!=0) amplitude[noos]*=-1;
	  amplitude[noos]*=-  /* i^2 */
	    2.0*(1.0/k+1.0/(state[nr][i]-k+plus)-1.0/state[nr][i]);
	  out_state[noos][0]=state[nr][0]+1;
	  out_state[noos][1]=k;
	  out_state[noos][2]=state[nr][i]-k+plus;
	  off[0]=1;
	  off[1]=i;
	  RRR(i-1,state[nr][0],nr);
	  for (j=1; j<=state[nr][0]-1; ++j) 
	    out_state[noos][j+2]=R[j-1];
	  ++noos;	    
	}
      if (state[nr][0]>2)
	{
	  amplitude[noos]=1;
	  if ((i+1)%2!=0) amplitude[noos]*=-1;
	  if ((state[nr][0]*i)%2!=0) amplitude[noos]*=-1;
	  amplitude[noos]*=-  /* i^2 */
	    2.0*(1.0/state[nr][mod(i-1+state[nr][0],state[nr][0])]
		-1.0/(state[nr][i]
		      +state[nr][mod(i-1+state[nr][0],state[nr][0])]+plus)
		 +1.0/state[nr][i]);
	  out_state[noos][0]=state[nr][0]-1;
	  out_state[noos][1]=state[nr][i]
	    +state[nr][mod(i-1+state[nr][0],state[nr][0])]+plus;
	  off[0]=2;
	  off[1]=i;
	  off[2]=mod(i+state[nr][0]-1,state[nr][0]);
	  RRR(i-1,state[nr][0],nr);
	  for (j=1; j<=state[nr][0]-2; ++j) 
	    out_state[noos][j+1]=R[j-1];
	  ++noos;	          
	}
    }
}


void create_states()
/* creates Fock states from K=<K_min>(even!) to K=<KK>+1 */ 
{
  int l;

  nos=0;
  /* -------------------- create fermion states */
  sector[K_min-3]=0;
  for (l=K_min-1; l<=KK+1; l+=2) 
    {
      find_state(1,l);
      sector[l]=nos;
    }
  nos2=nos;

  /* -------------------- create boson states */
  sector[K_min-2]=nos;
  for (l=K_min; l<=KK+2; l+=2) 
    {
      find_state(1,l);
      sector[l]=nos;
    }
}

void create_states_trunc(int p)
/* creates Fock states from K=<K_min>(even!) to K=<KK>+1, truncate to <p> */ 
/* partons */
{
  int l;

  nos=0;
  /* -------------------- create fermion states */
  sector[K_min-3]=0;
  for (l=K_min-1; l<=KK+1; l+=2) 
    {
      find_state_trunc(1,l,p+1);
      sector[l]=nos;
    }
  nos2=nos;

  /* -------------------- create boson states */
  sector[K_min-2]=nos;
  for (l=K_min; l<=KK+2; l+=2) 
    {
      find_state_trunc(1,l,p);
      sector[l]=nos;
    }
}


void construct_Q_minus(double *Q_Minus, int sign)
/* construct Q^-_{<sign>1/2} */
{
  int i,j,k,no,index=0;

  if (abs(sign)!=1) 
    {
      if (sign>0) sign=1;
      else sign=-1;
    }
  for (i=0; i<nos; ++i)           /* set Q_minus = 0.0 */
    for (j=0; j<nos; ++j) 
      {
	Q_Minus[index]=0.0;
	++index;
      }  
  for (k=0; k<nos; ++k)
    {
      noos=0;
      Q_minus(sign,k);
      if (test==1) printf("Q^-|%1d>=\n",k);
      for (i=0; i<noos; ++i)
	{
	  if (test==1)
	    {
	      printf("%8.5f |%1d; %1d",
		     amplitude[i],out_state[i][0],out_state[i][1]);
	      for (j=2; j<=out_state[i][0]; ++j) 
		printf(", %1d",out_state[i][j]);
	      printf(">\n");
	    }
	  no=get_state(i);
	  if (no>=0)
	    {
	      if ((perm_nr==0)||(out_state[i][0]%2!=0))
		Q_Minus[no+k*nos]+=amplitude[i]/norm[no]*norm[k];
	      else 
		{
		  if (perm_nr%2==0) 
		    Q_Minus[no+k*nos]+=amplitude[i]/norm[no]*norm[k];
		  else Q_Minus[no+k*nos]-=amplitude[i]/norm[no]*norm[k];
		}
	    }
	}
    }
}


void construct_Hamiltonian(double *P, double *Q1, double *Q2)
/* construct Hamiltonian P^-={Q_+1/2,Q_-1/2} */
{
  int i,j,k;

  printf("Hamiltonian:\n------------\n");
  for (i=0; i<nos; ++i)           /* construct Hamiltonian */
    {
      if ((test==1)&&(i==nos2)) printf("\n");
      for (j=0; j<nos; ++j) 
	{
	  if ((test==1)&&(j==nos2)) printf(" ");
	  /*P[i+j*nos]=0.0;*/	  
	  P[j+i*nos]=0.0;
	  for (k=0; k<nos; ++k)
	    /*P[i+j*nos]-=Q1[i+k*nos]*Q2[k+j*nos]+Q2[i+k*nos]*Q1[k+j*nos];*/ 
	    P[j+i*nos]-=Q1[k+i*nos]*Q2[j+k*nos]+Q2[k+i*nos]*Q1[j+k*nos]; 
	  if (test==1) printf("%5.3f ", P[i+j*nos]);
	}  
      if (test==1) printf("\n");
    }
  printf("\n");  
  printf("-----------------------------\n\n");
}


void construct_essential_Hamiltonian(double *ferm, double *bose)
/* construct Hamiltonian P^-={Q_+1/2,Q_-1/2} */
/* only relevant blocks: P^-=F^-_{K-1}*B^+_K + B^+_K*F^-_{K-1} */
{
  int i,j,k,dimF,dimB,dimB2,dimF2,no;
  long counter;
  double *F,*B,*F2,*B2,*F3,*B3;

  dimF=sector[KK-1];
  dimF2=nos2-sector[KK-1];
  dimB=sector[KK]-nos2;
  dimB2=sector[KK+2]-sector[KK];
  printf("Dimensions: 1st fermion block = %2d\n",dimF);
  F2=dvector(0,dimF2*dimB-1);
  printf("            2nd fermion block = %2d\n",dimF2);
  B=dvector(0,dimB*dimF-1);
  printf("            2nd boson block   = %2d\n",dimB);
  F=dvector(0,dimF*dimB-1);
  printf("            2nd boson block   = %2d\n",dimB2);
  B2=dvector(0,dimB*dimF2-1);
  if ((TRUNC<KK)&&(TruncType==1)) printf("Truncating B3 and F3...\n");
  else
    {
      F3=dvector(0,dimF2*dimB2-1);
      B3=dvector(0,dimB2*dimF2-1);
      for (i=0; i<dimF2; ++i)
	for (j=0; j<dimB2; ++j)
	  F3[i*dimB+j]=B3[i*dimB2+j]=0.0;
    }
  for (i=0; i<dimF; ++i)
    for (j=0; j<dimB; ++j)
      F[i*dimB+j]=B[i*dimB+j]=0.0;
  for (i=0; i<dimF2; ++i)
    for (j=0; j<dimB; ++j)
      F2[i*dimB+j]=B2[i*dimB+j]=0.0;
   
  for (k=nos2; k<sector[KK]; ++k)     /* construct F^-_{K-1} */
    {
      noos=0;
      Q_minus(-1,k);
      if (test==1) printf("Q^-|%1d>=\n",k);
      for (i=0; i<noos; ++i)
	{
	  if (test==1)
	    {
	      printf("%8.5f |%1d; %1d",
		     amplitude[i],out_state[i][0],out_state[i][1]);
	      for (j=2; j<=out_state[i][0]; ++j) 
		printf(", %1d",out_state[i][j]);
	      printf(">\n");
	    }
	  no=get_state(i);
	  if (no>=0)
	    {
	      if ((perm_nr==0)||(out_state[i][0]%2!=0))
		F[no+(k-nos2)*dimF]+=amplitude[i]/norm[no]*norm[k];
	      else 
		{
		  if (perm_nr%2==0) 
		    F[no+(k-nos2)*dimF]+=amplitude[i]/norm[no]*norm[k];
		  else F[no+(k-nos2)*dimF]-=amplitude[i]/norm[no]*norm[k];
		}
	    }
	}
    }

  for (k=nos2; k<sector[KK]; ++k)     /* construct F^+_{K+1} */
    {
      noos=0;
      Q_minus(1,k);
      if (test==1) printf("Q^-|%1d>=\n",k);
      for (i=0; i<noos; ++i)
	{
	  if (test==1)
	    {
	      printf("%8.5f |%1d; %1d",
		     amplitude[i],out_state[i][0],out_state[i][1]);
	      for (j=2; j<=out_state[i][0]; ++j) 
		printf(", %1d",out_state[i][j]);
	      printf(">\n");
	    }
	  no=get_state(i);
	  if (no>=0)
	    {
	      if ((perm_nr==0)||(out_state[i][0]%2!=0))
		F2[no-sector[KK-1]+(k-nos2)*dimF2]+=
		  amplitude[i]/norm[no]*norm[k];
	      else 
		{
		  if (perm_nr%2==0) 
		    F2[no-sector[KK-1]+(k-nos2)*dimF2]+=
		      amplitude[i]/norm[no]*norm[k];
		  else F2[no-sector[KK-1]+(k-nos2)*dimF2]-=
			 amplitude[i]/norm[no]*norm[k];
		}
	    }
	}
    }
  if ((TRUNC>=KK)||(TruncType==0))
    {
      for (k=sector[KK]; k<nos; ++k)     /* construct F^-_{K+1} */
	{
	  noos=0;
	  Q_minus(-1,k);
	  if (test==1) printf("Q^-|%1d>=\n",k);
	  for (i=0; i<noos; ++i)
	    {
	      if (test==1)
		{
		  printf("%8.5f |%1d; %1d",
			 amplitude[i],out_state[i][0],out_state[i][1]);
		  for (j=2; j<=out_state[i][0]; ++j) 
		    printf(", %1d",out_state[i][j]);
		  printf(">\n");
		}
	      no=get_state(i);
	      if (no>=0)
		{
		  if ((perm_nr==0)||(out_state[i][0]%2!=0))
		    F3[no-sector[KK-1]+(k-sector[KK])*dimF2]+=
		      amplitude[i]/norm[no]*norm[k];
		  else 
		    {
		      if (perm_nr%2==0) 
			F3[no-sector[KK-1]+(k-sector[KK])*dimF2]+=
			  amplitude[i]/norm[no]*norm[k];
		      else F3[no-sector[KK-1]+(k-sector[KK])*dimF2]-=
			     amplitude[i]/norm[no]*norm[k];
		    }
		}
	    }
	}
    }
  for (k=0; k<dimF; ++k)     /* construct B^+_{K-1} */
    {
      noos=0;
      Q_minus(1,k);
      if (test==1) printf("Q^-_+|%1d>=\n",k);
      for (i=0; i<noos; ++i)
	{
	  if (test==1)
	    {
	      printf("%8.5f |%1d; %1d",
		     amplitude[i],out_state[i][0],out_state[i][1]);
	      for (j=2; j<=out_state[i][0]; ++j) 
		printf(", %1d",out_state[i][j]);
	      printf(">\n");
	    }
	  no=get_state(i);
	  if (no>=0)
	    {
	      if ((perm_nr==0)||(out_state[i][0]%2!=0))
		B[no-nos2+k*dimB]+=amplitude[i]/norm[no]*norm[k];
	      else 
		{
		  if (perm_nr%2==0) 
		    B[no-nos2+k*dimB]+=amplitude[i]/norm[no]*norm[k];
		  else B[no-nos2+k*dimB]-=amplitude[i]/norm[no]*norm[k];
		}
	    }
	}
    }

  for (k=dimF; k<nos2; ++k)     /* construct B^-_{K} */
    {
      noos=0;
      Q_minus(-1,k);
      if (test==1) printf("Q^-_+|%1d>=\n",k);
      for (i=0; i<noos; ++i)
	{
	  if (test==1)
	    {
	      printf("%8.5f |%1d; %1d",
		     amplitude[i],out_state[i][0],out_state[i][1]);
	      for (j=2; j<=out_state[i][0]; ++j) 
		printf(", %1d",out_state[i][j]);
	      printf(">\n");
	    }
	  no=get_state(i);
	  if (no>=0)
	    {
	      if ((perm_nr==0)||(out_state[i][0]%2!=0))
		B2[no-nos2+(k-dimF)*dimB]+=amplitude[i]/norm[no]*norm[k];
	      else 
		{
		  if (perm_nr%2==0) 
		    B2[no-nos2+(k-dimF)*dimB]+=amplitude[i]/norm[no]*norm[k];
		  else B2[no-nos2+(k-dimF)*dimB]-=
			 amplitude[i]/norm[no]*norm[k];
		}
	    }
	}
    }

  if ((TRUNC>=KK)||(TruncType==0))
    {
      for (k=dimF; k<nos2; ++k)     /* construct B^+_{K+2} */
	{
	  noos=0;
	  Q_minus(1,k);
	  if (test==1) printf("Q^-_+|%1d>=\n",k);
	  for (i=0; i<noos; ++i)
	    {
	      if (test==1)
		{
		  printf("%8.5f |%1d; %1d",
			 amplitude[i],out_state[i][0],out_state[i][1]);
		  for (j=2; j<=out_state[i][0]; ++j) 
		    printf(", %1d",out_state[i][j]);
		  printf(">\n");
		}
	      no=get_state(i);
	      if (no>=0)
		{
		  if ((perm_nr==0)||(out_state[i][0]%2!=0))
		    B3[no-sector[KK]+(k-dimF)*dimB2]+=
		      amplitude[i]/norm[no]*norm[k];
		  else 
		    {
		      if (perm_nr%2==0) 
			B3[no-sector[KK]+(k-dimF)*dimB2]+=
			  amplitude[i]/norm[no]*norm[k];
		      else B3[no-sector[KK]+(k-dimF)*dimB2]-=
			     amplitude[i]/norm[no]*norm[k];
		    }
		}
	    }
	}
      counter=0;
      for (i=0; i<dimF2; ++i)   
	{
	  for (j=0; j<dimB2; ++j) 
	    {
	      if (fabs(B3[i*dimB2+j])>1e-8) ++counter;
	      /*printf("%5.3f ",B3[i*dimB2+j]);*/ 
	    }
	  /*printf("\n");*/
	}
      printf("Non-zero in B3: %2d = %5.3f percent\n",
	     counter,1.0*counter/dimF2/dimB2);


      for (i=0; i<dimF2; ++i)       /* construct fermionic Hamiltonian */
	{
	  if ((test==1)&&(i==nos2)) printf("\n");
	  for (j=0; j<dimF2; ++j) 
	    {
	      if ((test==1)&&(j==nos2)) printf(" ");
	      ferm[i+j*dimF2]=0.0;
	      for (k=0; k<dimB2; ++k) 
		ferm[i+j*dimF2]-=B3[k+i*dimB2]*F3[j+k*dimF2];
	      for (k=0; k<dimB; ++k) 
		ferm[i+j*dimF2]-=B2[k+i*dimB]*F2[j+k*dimF2];
	      if (test==1) printf("%5.3f ", ferm[i+j*dimF]);
	    }  
	  if (test==1) printf("\n");
	}
    }
  printf("Constructing bosonic part of Hamiltonian...\n");  
  for (i=0; i<dimB; ++i)       /* construct bosonic Hamiltonian */
    {
      if ((test==1)&&(i==nos2)) printf("\n");
      for (j=0; j<dimB; ++j) 
	{
	  if ((test==1)&&(j==nos2)) printf(" ");
	  bose[i+j*dimB]=0.0;
	  for (k=0; k<dimF; ++k) 
	    bose[i+j*dimB]-=F[k+i*dimF]*B[j+k*dimB];
	  for (k=0; k<dimF2; ++k) 
	    bose[i+j*dimB]-=F2[k+i*dimF2]*B2[j+k*dimB];
	  if (test==1) printf("%5.3f ", bose[i+j*dimB]);
	}  
      if (test==1) printf("\n");
    }
  printf("\n");  
}

void get_intact_Hamiltonian(double *fmat, double *bmat, int n)
/* chop off of Hamiltonian <Ham> the incomplete parts, */
/* i.e. the first and last block of both fermions and bosons, */
/* puts result in <fmat> and <bmat> */
{
  int i,j,k,index=0;

  for (k=K_min-1; k<=KK+1; k+=2) /* loop over all sectors */
    {
      for (i=sector[k-2]; i<sector[k]; ++i)
	{
	  for (j=0; j<sector[k-2]; ++j) ++index;
	  for (j=sector[k-2]; j<n; ++j)
	    {
	      if ((j<sector[k])&&(i>=sector[K_min-1])) 
		{
		  fmat[(i-sector[K_min-1])*(sector[KK+1]-sector[K_min-1])
		      +j-sector[K_min-1]]=Ham[index];
		}
	      ++index;
	    }
	}
    }

  /* bosons */
 
  for (k=K_min; k<=KK; k+=2) /* loop over all sectors except last one */
    {
      for (i=sector[k-2]; i<sector[k]; ++i)
	{
	  for (j=0; j<sector[k-2]; ++j) ++index;
	  for (j=sector[k-2]; j<n; ++j)
	    {
	      if ((j<sector[k])) 
		{
		  bmat[(i-sector[KK+1])*(sector[KK]-sector[KK+1])+j
		      -sector[KK+1]]=Ham[index];
		}
	      ++index;
	    }
	}
    }
}


void print_eigenvalues(int dim, int fb)
/* prints eigenvalues, dimension <dim> */ 
/* <fb>=0 --> boson */
/* output file "evs_<fb>K<KK>.dat" */
{
  int i,shift=0;
  char *name;
  FILE *fp;

  name=(char *)malloc(80);
  if (fb==1) shift=1;
  sprintf(name,"evs_%1dK%1d.dat",fb,KK+shift);
  fp=fopen(name,"w");

  cout<<"Writing eigenvalues into <"<<name<<">"<<endl;

  for (i=0; i<dim; ++i) fprintf(fp,"%2d  %12.8f\n",KK+shift,r[i]);    
  
  fclose(fp);
}

void print_eigenfunctions(int nr, double *mat, int dim, int fb)
/* prints eigenfunctions up to <nr>-1 stored in <mat> of dimension <dim> */ 
/* <fb>=0 --> boson */
/* output file "efs_<fb>K<KK>.dat" */
{
  int i,j,k,count=0,shift;
  char *name;
  FILE *fp;

  name=(char *)malloc(80);
  if (fb==1) sprintf(name,"efs_%1dK%1d.dat",fb,KK+1);
  else sprintf(name,"efs_%1dK%1d.dat",fb,KK);
  fp=fopen(name,"w");
  cout<<"Writing eigenfunctions into <"<<name<<">"<<endl;

  if (fb==1) 
    {
      shift=sector[K_min-1];
      fprintf(fp,"Fermionic eigenfunctions at K=%2d:\n\n",KK+1);
    }
  else 
    {
      shift=sector[KK+1];
      fprintf(fp,"Bosonic eigenfunctions at K=%2d:\n\n",KK);
    }
  for (k=0; k<nr; ++k)
    {
      fprintf(fp,"|%1d.ef, M^2=%5.3f>=\n",k,r[k]);
      for (i=0; i<dim; ++i)
	{
	  if (fabs(mat[i*dim+numb[k]])>ef_relevance)
	    {
	      if (count%3==0) fprintf(fp,"\n");
	      fprintf(fp,"+(%5.3f)|%1d; %1d",
		   mat[i*dim+numb[k]],state[i+shift][0],state[i+shift][1]);
	      for (j=2; j<=state[i+shift][0]; ++j) 
		fprintf(fp,",%1d",state[i+shift][j]);
	      fprintf(fp,">");
	      ++count;
	    }
	}  
      fprintf(fp,"\n\n");
    }
  fclose(fp);  
}

void sort_evs(int dim, double *evs, int *nr)
/* sort eigenvalues <nr> of dimension <dim>, store eigenfunction nr. in <nr>*/ 
{
  int i,j,count=0,dim2,lowj;
  double low;

  dim2=0;
  for (i=0; i<dim; i++) nr[i]=i;
  while (dim2<dim)
    {
      low=evs[dim2];
      lowj=nr[dim2];
      j=dim2;      
      for (i=dim2; i<dim; ++i)
	{
	  if (evs[i]<low) 
	    {
	      low=evs[i];
	      lowj=nr[i];
	      j=i;
	    }
	}
      for (i=j; i>dim2; --i) 
	{
	  evs[i]=evs[i-1];
	  nr[i]=nr[i-1];
	}
      evs[dim2]=low;
      nr[dim2]=lowj;
      ++dim2;
      low=evs[dim2];
    }
}


void main()
{
  int i,ii,j,fdim,bdim;
  double *fermions,*bosons;

  cout<<endl<<" S Y M _ w i t h _ a p B C . c p p      Version Aug 2001"<<endl;
  cout<<"==========================================================="<<endl;
  cout<<endl<<"   Calculates the spectrum of SYM(1+1) with"<<endl;
  cout<<"   anti-periodic boundary conditions, see paper "<<endl;      
  cout<<"   S.Pinsky, U.Trittmann, Phys.Rev. D62 (2000) 087701"<<endl;
  cout<<"   (hep-th/0005055)."<<endl<<endl;
  cout<<"-----------------------------------------------------------"<<endl;
  cout<<endl<<" K_max     = "<<KKK<<endl;
  cout<<" TRUNC     = "<<TRUNC<<endl;
  cout<<" TruncType = "<<TruncType<<endl;
  cout<<endl;


  for (KK=4; KK<=KKK; KK+=2)
    {
      printf("\nK/K_max = %2d/%2d \n-----------------------------\n",KK,KKK);
      test=0;
      K_min=KK;
      create_states_trunc(TRUNC);
      test=0;
      
      fdim=sector[KK+1]-sector[KK-1];
      bdim=sector[KK]-nos2;
      if ((TRUNC>=KK)||(TruncType==0)) fermions=dvector(0,fdim*fdim-1);
      else fermions=dvector(0,1);
      bosons=dvector(0,bdim*bdim-1);
      
      test=0;
      printf("Constructing Hamiltonian...\n");
      construct_essential_Hamiltonian(fermions,bosons);
      test=0;
      
      if ((TRUNC>=KK)||(TruncType==0)) 
	{
	  diagonalize_symmetric(&fermions,fdim,&r,&v);
	  for (i=0; i<fdim; ++i) r[i]*=(KK+1)/4.0;
	  numb=intvector(0,fdim);
	  sort_evs(fdim,r,numb);

	  ii=6;
	  if (fdim<6) ii=fdim;      
	  printf("\n Fermions:\n");
	  for (i=0; i<ii; ++i) printf("%2d %12.8f\n",i,r[i]);
      
	  print_eigenfunctions(ii,fermions,fdim,1);
	  print_eigenvalues(fdim,1);
	  free_Intvector(numb,0,fdim-1);
      
	  free_dvector(r,0,fdim-1);
	}

      diagonalize_symmetric(&bosons,bdim,&r,&v);
      for (i=0; i<bdim; ++i) r[i]*=KK/4.0;
      numb=intvector(0,bdim);
      sort_evs(bdim,r,numb);
      
      ii=6;
      if (bdim<6) ii=bdim;      
      printf("\nBosons:\n");
      for (i=0; i<ii; ++i) printf("%2d %12.8f\n",i,r[i]);
      
      print_eigenfunctions(ii,bosons,bdim,0);
      print_eigenvalues(bdim,0);
    }
}