#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "basis.h"

/* 
         Subroutines for generating bases of states. 
         Summer 2002 version, including mesons.

         if state[...] is an array holding a meson state.
         state[0] is the quark and state[np-1] is the antiquark.

	 The HEAVY_LIGHT states consist of a heavy quark
         and a light anti-quark.  This is the assumption used
         for determining any transformation properties.
 
         The SOURCE states assume that the heavy fields
         transform as scalars.
*/

#define PRINTIT 0 /*print basis size & timing information in makebasis */

/*  initbasis(x,type,ht,nptrunc,kt,charge,o,multi) 
    initializes the basis parameters in x.

   If o == 0, no orientation reversal symmetry
      o == 1, even symmetry
      o ==-1, odd symmetry
    Likewise for charge conjugation.

    kt is twice the total discretized momentum.
*/

void initbasis(full_basis *x, enum state_type type, int ht[HINDEX], 
	       int nptrunc, unsigned int kt, int charge, int o, int multi){
  int i;
  if(x->length>0)free_basis(x);
  x->type=type;
  for(i=0; i<HINDEX; i++)x->ht[i]=ht[i];
  x->nptrunc=nptrunc; x->kt=kt;  
  x->o=o; x->charge=charge; x->multi=multi;
}

/* 
   Free storage allocated for full_basis x.
*/

void free_basis(full_basis *x){
  free(x->inner);
  free(x->np);
  free(x->s);
  x->length=0;
}

/* makebasis(x) finds the inner product list and the actual
     basis of states.  */

void makebasis(full_basis *x){
  int np,i,ss;
  size_t ist;
  partype *pointx,*pointy;
  reduced_basis *y;
#if PRINTIT
  clock_t t1=clock();
#endif

  /* Complain if a basis has already been constructed. */
  if(x->length>0){
     fprintf(stderr,__FILE__":  Basis already filled in makebasis\n");
    goto error;
  }  

  makereduced(&y,x->type,x->ht,x->nptrunc,x->kt,x->charge,x->o,
	      subgroup[x->multi]);

/*
  Set aside memory for a copy of the basis along with
  other information (like the inner product).
*/

  pointx=x->s=malloc(y->length*x->nptrunc*sizeof(*(x->s)));
  x->np=malloc(y->length*sizeof(*(x->np)));
  x->inner=malloc(y->length*sizeof(*(x->inner)));

  for(ist=0, pointy=y->s; ist<y->length; ist++, pointy+=y->nptrunc){
    for(np=0; np<y->nptrunc && pointy[np]!=HOLE; np++);
    ss=innerproduct(pointy,pointy,y->type,np,x->charge,x->o,x->multi);
    if(ss){
      memcpy(pointx,pointy,x->nptrunc*sizeof(partype));
      pointx+=x->nptrunc;
      x->np[x->length]=np;
      x->inner[x->length]=ss;
      x->length++;
    }
  }

  if(x->length==0){
     fprintf(stderr,__FILE__":  Empty basis\n");
    goto error;
  }

  /*  Print out the results */

#if 0  /* debugging stuff */
  printf("    The basis is (with inner product):\n");
  for(ist=0; ist<x->length; ist++){
    printf("     ");
    printstate((x->s)+ist*(x->nptrunc),x->type,x->nptrunc);
    printf("  %i\n",x->inner[ist]);
  }
#endif
#if PRINTIT
  printf("      Basis: %i nonzero norm states",(int) x->length);
  printf(" in %.2f CPU sec\n",
	 (float)(clock()-t1)/(float) CLOCKS_PER_SEC);
#endif
  return;

error:
  fprintf(stderr,__FILE__":  Error in makebasis, ht=");
  for(i=0; i<HINDEX; i++) fprintf(stderr," %i",x->ht[i]);
   fprintf(stderr,"  , np=%i, K=%i/2, charge=%i, o=%i, multi=%i, exiting\n",
	  x->nptrunc,x->kt, x->charge, x->o, x->multi);
  free_basis(x);
  return;
}



/*  give state canonical order with flags for
    orientation reversals and charge conjugation symmetries,
    and transverse lattice symmetry subgroup
    (see group.c for definitions).

    It is not a good idea to use memcmp() in the following 
    because if partype spans two bytes, memcmp() can assign
    an incorrect order for the two states.

*/

#define COMPARE(Z)                                       \
  for(i=0; i<np; i++){                                   \
    if(best[i]==(Z)[i])continue;                         \
    if(best[i]>(Z)[i])memcpy(best,Z,np*sizeof(partype)); \
    break;                                               \
  }                                                      \

void wellorder(partype *state, enum state_type type, int np, 
	       int charge, int o, int subgroup){
  int i,cycle;
  partype best[NPMAX],z1[2*NPMAX],z2[2*NPMAX];
  int g;   
  const latlist *subg=latgroups+subgroup;
  memcpy(best,state,np*sizeof(partype));
  for(g=0; g<subg->length; g++){
    statelat(z1,state,type,np,subg->e[g]);
    /* can't have both "charge" and "o" */
    if(charge)
      charge_conjugation(z2,z1,type,np);
    else if(o)
      orientation(z2,z1,type,np);
    /*  This avoids the use of the "%" which is really slow */
    /*  For speed, moved function calls out of the cycle loop */
    switch(type){
    case TYPE_GLUE:
      memcpy(z1+np,z1,np*sizeof(partype));
      memcpy(z2+np,z2,np*sizeof(partype));
      for(cycle=0; cycle<np; cycle++){
	/* pointer for cyclic permutations */
	COMPARE(z1+cycle);
	if(charge||o)COMPARE(z2+cycle);
      }
      break;
    default:
      COMPARE(z1);
      if(charge||o)COMPARE(z2);
      break;
    }
  }
  memcpy(state,best,np*sizeof(partype));
}

/* 
    makereduced makes a basis of states

    If rotates nonzero helicities to a canonical direction.  This
    will give a factor of HINDEX! * 2^HINDEX improvement in memory and 
    speed. 

*/ 

#define USE_TREE 1 /* whether to use binary tree to construct basis */

#define MAX_BASE 10000 /* number of possible bases */
static reduced_basis basislist[MAX_BASE]; /* Table of bases */
static int basisfilled=0;                 /* Number of bases in basislist */
size_mem reduced_basis_memory=0;  /* tally memory used to store bases */

int makereduced(reduced_basis **x, enum state_type type, int ht[HINDEX], 
		int nptrunc, unsigned int kt, int charge, int o, int subgroup){

  int i,j,l,hlist[HINDEX],abswind=0,reverse, r_group, cycle;
  int kk,hh1,shift;  /* these must be signed variables */
  int two_trip;
  int temp;
  int nptrunc_less;
  int hlink[HINDEX];
  int g;
  latype bestdirection=0; /* Lattice transform for gluon states */
  partype hh;
  size_t length=0;
  partype *s=NULL,z[2*NPMAX],z2[NPMAX];
  partype *pointy, particle;
  reduced_basis *y;
  int npmin=999; 
  const latlist *subg=latgroups+subgroup;
#if USE_TREE
  tree_basis *tree=NULL;
#endif

/*
       Check if it has already been done 
*/

  for(y=basislist, i=0; i<basisfilled; y++, i++){
    if(y->type != type) continue;
    if(type==TYPE_GLUE || type==TYPE_SOURCE){
      for(j=0; j<HINDEX && y->ht[j]==ht[j]; j++);
      if(j<HINDEX)continue;
    }
    if(y->kt==kt && y->charge==charge && y->o==o &&
       subgroup == y->subgroup && y->nptrunc==nptrunc){
      *x=y;
#if 0 /* debugging print */
       fprintf(stderr,__FILE__":    Existing: ht =");
      for(i=0;i<HINDEX;i++) fprintf(stderr," %2i",ht[i]);
       fprintf(stderr,", nptrunc = %i, K = %i/2\n", nptrunc, kt); 
#endif 
      return 0;
    }
  }

/* 
   Minimum number of particles allowed in state
*/
  switch(type){
  case TYPE_MESON: 
    npmin=2; 
    break;
  case TYPE_HEAVY_LIGHT:
  case TYPE_GLUE:
    npmin=1; 
    break;
  case TYPE_SOURCE:
    npmin=0;
    break;
  default:
    fprintf(stderr,__FILE__":  invalid state type\n");
    return -1;
  }
  
  /*   
       Check various error conditions.
  */
  
#if 0 /* debugging print */
   fprintf(stderr,__FILE__":    New:  ht =");
  for(i=0;i<HINDEX;i++) fprintf(stderr," %2i",ht[i]);
   fprintf(stderr,", nptrunc = %i, K = %i/2\n", nptrunc, kt); 
#endif 
  
  if(basisfilled==MAX_BASE){
    fprintf(stderr,__FILE__":  MAX_BASE too small, exiting\n");
    return -1;
  } 
  
  if(nptrunc>NPMAX || nptrunc<npmin){
    fprintf(stderr,__FILE__":  nptrunc=%i is not between %i and NPMAX=%i\n",
	    nptrunc,npmin,NPMAX);
    return -1;
  }
  
  if(charge && o){
    fprintf(stderr,__FILE__":  specifying orientation=%i and charge" 
	    " conjugation=%i is redundant\n",o,charge);
    return -1;
  }
  
  /* check that the winding number is appropriate */
  if(type==TYPE_GLUE || type==TYPE_SOURCE){  
    abswind=0;
    for(l=0; l<HINDEX; l++)abswind += abs(ht[l]);
    if(abswind>nptrunc || nptrunc%2!=abswind%2){
      fprintf(stderr,__FILE__":  illegal value for ht[], nptrunc=%i\n",
	      nptrunc);
      return -2;
    }
    if(abswind && charge){
      fprintf(stderr,__FILE__":  can't have charge conjugation symmetry\n"
	      "               and nonzero winding, abswind=%i\n",abswind);
      return -3;
    }
  }
  
#if !GLUE_PERIODIC
  if(type==TYPE_MESON || type==TYPE_HEAVY_LIGHT){
    fprintf(stderr,__FILE__":  no mesons for antiperiodic link fields\n");
    return -1;
  }
#endif
  
  /*
    Make a basis, first the few-particle case
  */
  
  if(nptrunc==npmin)
    switch(type){
    case TYPE_GLUE:
      if(link_k(kt)){
	length=1;
	s=realloc(s,length*nptrunc*sizeof(partype));
	s[0]=glueencode(kt/2,ht);
      }
      break;
      
    case TYPE_MESON:
      s=realloc(s,4*kt/2*nptrunc*sizeof(partype)); 
#if !QUARK_PERIODIC /* didn't write a periodic version */
      if(kt%2==1)break;  /* improper value for momentum */
      for(kk=0; kk<kt/2; kk++) /* kk is the quark momentum rounded down */
	for(i=0; i<2; i++) 
	  for(j=0; j<2; j++){ 
	    z[0]=quarkencode(i,kk);
	    z[1]=quarkencode(j,kt/2-kk-1);
#if 0  /* debug print */
	    printf(" (%i,%i) (%i,%i)\n",z[0]>>HSHIFT,z[0]&KMASK,
		   z[1]>>HSHIFT,z[1]&KMASK);
#endif
	    wellorder(z,type,nptrunc,charge,o,subgroup);
	    length=addstate(s,z,length,nptrunc);
	  }
#else
      fprintf(stderr,__FILE__":  Periodic for quarks\n");
      return -1;
#endif
      s=realloc(s,length*nptrunc*sizeof(partype));
      break;
      
    case TYPE_HEAVY_LIGHT:
      s=realloc(s,(kt+1)*nptrunc*sizeof(partype));
#if !QUARK_PERIODIC /* didn't write a periodic version */
      for(kk=0; 2*kk+1<=kt; kk++) /* kk is the quark momentum rounded down */
	for(i=0; i<2; i++){ 
	  z[0]=quarkencode(i,kk);
#if 0  /* debug print */
	  printf(" (%i,%i) (%i,%i)\n",z[0]>>HSHIFT,z[0]&KMASK,
		 z[1]>>HSHIFT,z[1]&KMASK);
#endif
	  wellorder(z,type,nptrunc,charge,o,subgroup);
	  length=addstate(s,z,length,nptrunc);
	}
#else
      fprintf(stderr,__FILE__":  Periodic for quarks\n");
      return -1;
#endif
      s=realloc(s,length*nptrunc*sizeof(partype));
      break;
      
    case TYPE_SOURCE:
      length=1;
      break;
    } 
  
  else {  /******   Many particle case *************/
    
    /*  Fill basis with states containing fewer particles */
    /* abswind defined above for non-meson states */
    if(type==TYPE_MESON || type==TYPE_HEAVY_LIGHT)
      nptrunc_less=nptrunc-1;
    else
      nptrunc_less=nptrunc-2;
    if(nptrunc_less>=npmin && 
       (type==TYPE_MESON || type==TYPE_HEAVY_LIGHT || nptrunc_less>=abswind)){ 
      makereduced(&y,type,ht,nptrunc_less,kt,charge,o,subgroup);
#if USE_TREE
      tree=tree_realloc(tree,length+y->length,nptrunc);
#else
      s=realloc(s,(length+y->length)*nptrunc*sizeof(partype));
#endif
#if 1  /* for binary tree, go through basis out of order */
      for(shift=1; y->length>=shift; shift*=2);
      while(shift>>=1 >= 1)
	for(j=shift-1; j<y->length; j+=2*shift)
#else  /* go through basis in order */
	for(j=0; j<y->length; j++)
#endif
	{
	  pointy=y->s+y->nptrunc*j;
	  for(i=0; i<y->nptrunc; i++)
	    z[i]=pointy[i];
	  for(i=y->nptrunc; i<nptrunc; i++)
	    z[i]=HOLE;  /* fill empty slots with HOLE */
	  /* below, the last empty slot must have HOLE */
#if USE_TREE
	  length=tree_addstate(tree,z,length,nptrunc);
#else
	  length=addstate(s,z,length,nptrunc);
#endif
	}
    }
    
    /*    
	  nptrunc particle case: add a link field to the state 
	  
	  The kk and hh1 loops are done backwards to increase
	  the randomization of the order of the states.
	  This is needed for tree_addstate to work effectively.
	  
	  in order to reduce "holes" in memory, we go through
	  the loops twice.  The first time, we find the bases with
	  fewer particles, and the second time we actually calculate
	  the states in the new basis.
    */
    
    for(two_trip=0; two_trip<2; two_trip++)    
#if 1 /* done backwards to increase speed for binary tree */
      for(hh1=2*HINDEX-1; hh1>=0; hh1--)
#else
      for(hh1=0; hh1<2*HINDEX; hh1++)
#endif
      {
	/* hh can not be used as the loop variable because
	   of possible overflow when hh = (2*HINDEX)<<HSHIFT  */
	hh=hh1<<HSHIFT;
	(void) gluedecode(hh,hlink);
	
	/*  Make sure that we haven't already used this index assuming
	    equivalence under the specified subgroup.            */
	for(g=0; g<subg->length; g++){
	  partype hhtemp;
	  statelat(&hhtemp,&hh,TYPE_SOURCE,1,subg->e[g]);
	  if(hh>hhtemp)break;
	}
	/* if we found an hh done previously, go to next hh. */
	if(g<subg->length)continue;
	
	for(kk=0; kk<=kt; kk++){  /* kk is twice the discretized momentum */
	  if(!link_k(kk)) continue; /* skip illegal values of momentum */
	  
	  switch(type){
	  case TYPE_MESON:
	    makereduced(&y,type,NULL,nptrunc-1,kt-kk,1,0,r_group=FULL_GROUP);
	    break;
	  case TYPE_HEAVY_LIGHT:
	    makereduced(&y,type,NULL,nptrunc-1,kt-kk,0,0,r_group=FULL_GROUP);
	    break;
	  default:
	    /* Find the minimum number of particles necessary 
	       for the winding number */
	    {int abswind_temp=0;
	    for(l=0; l<HINDEX; l++)
	      abswind_temp += abs(hlist[l]=ht[l]-hlink[l]);
	    if(abswind_temp > nptrunc-1) continue;}
	    
	    /*  Change winding to a preferred direction */
	    bestdirection=latorder(hlist);
	    /*
	      we assume no lattice symmetries for
	      states with fewer particles since they have, in general, 
	      nonzero winding.  
	      However, we assume orientation reversal symmetry.  
	    */
	    
	    makereduced(&y,type,hlist,nptrunc-1,kt-kk,0,1,
			r_group=TRIVIAL_GROUP);
#if 0  /* debugging print */
	    printf("bestdirection=%i ",bestdirection);
	    printf("ht=(%i,%i) hlink=(%i,%i) hlist=(%i,%i)\n",ht[0],ht[1],hlink[0],
		   hlink[1],hlist[0],hlist[1]);
#endif
	  }
	  
	  if(two_trip==0)continue; /* Don't do anything the first time */
	  
	  /*  
	      add a particle to each state in the basis.
	      Undo all symmetries in basis states:
	      any lattice symmetries (r_group), orientation reversals,
	      or cyclic symmetries (for glueballs)
	  */	  
	  
	  particle=glueencode(kk/2,hlink);
	  
#if 1  /* for binary tree, go through basis out of order */
	  for(shift=1; y->length>=shift; shift*=2);
	  while(shift>>=1 >= 1)
	    for(i=shift-1; i<y->length; i+=2*shift)
#else  /* go through basis in order */
	    for(i=0; i<y->length; i++)
#endif
	    {
	      pointy=y->s+y->nptrunc*i;
	      /* check that the state actually has y->nptrunc particles */
	      if(y->nptrunc>0 && pointy[y->nptrunc-1]==HOLE)continue;
#if 1     /* For glueballs, only insert a link with equal or greater K */
	      /* this greatly increases the speed */
	      if(type==TYPE_GLUE){
		for(l=0; l<y->nptrunc; l++)
		  if((particle&KMASK)>(pointy[l]&KMASK))break;
		if(l<y->nptrunc)continue;
	      }
#endif
	      /* make sure there is enough memory available */
	      temp=length+y->length*latgroups[r_group].length*2*
		(type==TYPE_GLUE?y->nptrunc:1);
#if USE_TREE
	      tree=tree_realloc(tree,temp,nptrunc);
#else
	      s=realloc(s,temp*nptrunc*sizeof(partype));   
	      if(s==NULL){
		fprintf(stderr,__FILE__
			":  2 no memory left in makereduced, exiting\n");
		return -8;
	      } 
#endif
	      for(g=0; g<latgroups[r_group].length; g++){
		/* Generate basis states without lattice symmetries */
		statelat(z,pointy,type,y->nptrunc,latgroups[r_group].e[g]);
		/* For nonzero winding, change direction back */
		if(type==TYPE_GLUE || type==TYPE_SOURCE)
		  statelat(z,z,type,y->nptrunc,bestdirection);
		reverse=1; 
		do{
		  if(type==TYPE_GLUE){
		    memcpy(z+y->nptrunc,z,y->nptrunc*sizeof(partype));
		    cycle=y->nptrunc-1;
		  } else
		    cycle=0;
		  do{
		    for(l=0; l<y->nptrunc; l++)                  
		      z2[l]=z[l+cycle];                              
		    if(type==TYPE_MESON || type==TYPE_HEAVY_LIGHT){           
		      z2[nptrunc-1]=z[y->nptrunc-1];           
		      z2[nptrunc-2]=particle;                    
		    } else                                       
		      z2[nptrunc-1]=particle;          
#if GLUE_PERIODIC                               
		    if(!zeromodetest(z2,type,nptrunc))continue;  
#endif                                          
		    wellorder(z2,type,nptrunc,charge,o,subgroup);
#if USE_TREE
		    length=tree_addstate(tree,z2,length,nptrunc);
#else
		    length=addstate(s,z2,length,nptrunc);        
#endif
		  } while(cycle--);
		  /* undo any orientation reversal/charge conjugation symmetry
		     of fewer particle basis */ 
		  if(reverse){
		    if(type==TYPE_HEAVY_LIGHT)
		      reverse=0;
		    else if(type==TYPE_MESON)
		      charge_conjugation(z,z,type,y->nptrunc);
		    else
		      orientation(z,z,type,y->nptrunc);
		  }
		} while(reverse--);
	      }
	    }
    	}
      }
#if 0 /* debugging print */ 
    fprintf(stderr,"length = %i K=%i/2 nptrunc=%i\n",length,kt,nptrunc);
#endif
#if USE_TREE
    s=tree_sort(tree,length,nptrunc);
#else
    if(length*nptrunc>0)
      s=realloc(s,length*nptrunc*sizeof(partype));
    /* ANSI standard does not guarantee that realloc sets s to NULL */
    else {  
      free(s);
      s=NULL;
    }
#endif
  }

  /* tally memory used to store bases, in bytes */
  reduced_basis_memory += nptrunc*length*sizeof(partype);
  
  /* Here we store the results */
  *x=basislist+(basisfilled++);
  (*x)->type = type;
  if(type==TYPE_MESON || type==TYPE_HEAVY_LIGHT)
    (*x)->ht[0]=999;
  else
    for(i=0; i<HINDEX; i++) (*x)->ht[i]=ht[i];
  (*x)->kt = kt;
  (*x)->charge = charge;
  (*x)->o = o;
  (*x)->subgroup = subgroup;
  (*x)->nptrunc = nptrunc;
  (*x)->length = length;
  (*x)->s = s;
  return 0;
}

/* 
   allowed values of link momentum (true if allowed)
   kt is twice the discretized momentum.
*/

int link_k(unsigned int kt){  
  return 
    /* no illegal zero modes */
#if GLUE_PERIODIC
    !(kt==0 && ADJACENT==0) &&
#endif
    /* kt is even for periodic and odd for antiperiodic */
    kt%2==(GLUE_PERIODIC?0:1) &&
    /* Don't accidentally make a HOLE */  
    (kt>>1) < (1<<HSHIFT)-1;
}

/* 
    Determine whether a state has the allowed number adjacent
    link field zero modes.  Returns true if state is allowed.
*/

int zeromodetest(partype *s, enum state_type type, int np){

#if GLUE_PERIODIC
  int i,j=0;
 switch(type){
  case TYPE_GLUE:
    for(i=0; j<np && (i<np || j!=0); i++){
      if((s[i%np]&KMASK) == 0)
	j++;
      else 
	j=0;
      if(j>ADJACENT)return 0;
    }
    break;
  case TYPE_MESON:
    for(i=1; i<np-1; i++){
      if((s[i]&KMASK) == 0)
	j++;
      else 
	j=0;
      if(j>ADJACENT)return 0;
    }
    break;
  case TYPE_HEAVY_LIGHT:
    for(i=0; i<np-1; i++){
      if((s[i]&KMASK) == 0)
	j++;
      else 
	j=0;
      if(j>ADJACENT)return 0;
    }
    break;
  case TYPE_SOURCE:
    for(i=0; i<np; i++){
      if((s[i]&KMASK) == 0)
	j++;
      else 
	j=0;
      if(j>ADJACENT)return 0;
    }
    break;
  }
#endif
  return 1;
}

/*  
      Free up memory used by makereduced.
*/

void free_reduced(void){
  while(basisfilled>0)
      free(basislist[--basisfilled].s);
  reduced_basis_memory=0;
}

/* return a flag to say if a basis is winding  */

int iswinding(int *ht){
  int i;
  for(i=0; i<HINDEX && ht[i]==0; i++);
  return i<HINDEX;
}

/*
    multiply complex number by i^(phase/2)
*/

void phasemultiply(complex_element *z,int phase){
  complex_element temp;
  if(phase&1){
    temp.re=z->re;
    temp.im=z->im;
    z->re=(temp.re-temp.im)/sqrt(2.0);
    z->im=(temp.im+temp.re)/sqrt(2.0);
  }
  if(phase&2){
    temp.im=z->im;
    z->im=z->re;
    z->re=-temp.im;
  }
  if(phase&4){
    z->re*=-1;
    z->im*=-1;
  }
}

/*  Find the inner product of two states, with  
   orientation reversals, charge conjugation  and
   transverse lattice symmetry multiplet multi 
   (see group.c for definitions).  
   If o == 0, no orientation reversal symmetry
      o == 1, even symmetry
      o ==-1, odd symmetry
   Likewise for charge conjugation.
   This returns the number of nonzero contractions weighted by symmetries.

   This assumes the phase factor is -1, i, -i, or +1
	 There is an implicit mod 8
*/

#define INNER_COMPARE(Z,PHASE)  {                        \
  for(i=0; i<np; i++)                                    \
    if(left[i]!=(Z)[i])break;                            \
  if(i==np){                                             \
    if(PHASE&1) goto error;                              \
    else if(PHASE&3) s.im+=(PHASE&4?-1:1);               \
    else s.re+=(PHASE&4?-1:1);                           \
  }                                                      \
}                                                        \

int innerproduct(partype *left, partype *right, enum state_type type, int np, 
		 int charge, int o, int multi){
  struct {int re; int im;} s={0,0};  /* accumulator for weights */
  int i,cycle,phase=0,phase_reverse=999;
  partype z1[2*NPMAX],z2[2*NPMAX];
  int g;
  const latlist *subg=latgroups+subgroup[multi];
  
  for(g=0; g<subg->length; g++){
    phase=(multiplet[multi][g]==-1?4:0);
    phase+=statelat(z1,right,type,np,subg->e[g]);
    if(charge){
      phase_reverse=phase+charge_conjugation(z2,z1,type,np);
      if(charge==-1)phase_reverse+=4;
    } else if(o){
      phase_reverse=phase+orientation(z2,z1,type,np);
      if(o==-1)phase_reverse+=4;
    }
    
    if(type==TYPE_GLUE){
      memcpy(z1+np,z1,np*sizeof(partype));
      memcpy(z2+np,z2,np*sizeof(partype));
      for(cycle=0; cycle<np; cycle++){
	INNER_COMPARE(z1+cycle,phase);
	if(charge||o)INNER_COMPARE(z2+cycle,phase_reverse);
      }
    } else {
      INNER_COMPARE(z1,phase);
      if(charge||o)INNER_COMPARE(z2,phase_reverse);
    }
    
#if 0  /* debug flag */
    printf("      ");
    printstate(z,type,np);
    printf(" phase=i^%i\n",(phase%8)>>1);
#endif
  }
  if(s.im)goto error; /* inner is real */
  return s.re;

 error:
  fprintf(stderr,__FILE__
	  ":  Complex number in inner product! phase=%i\n",phase);
  fprintf(stderr,"  np=%i, charge=%i, o=%i, multi=%i\n  ",
	  np,charge,o,multi);
  for(i=0; i<np; i++)
    fprintf(stderr,"(%i,%i) ",left[i]>>HSHIFT,left[i]&KMASK);
  fprintf(stderr,"\n");
  exit(-1);
}


/*  
     lattice group element acting on a quark or anti-quark 
     creation operator.
     The phase factor associated with the operation is
     equal to  i^(phase/2).
     Note that the function *adds* to the phase.
*/

#if HINDEX==2
const int flip[ORDER_GROUP]={0,1,1,0,1,0,0,1};
#else
const int *flip=NULL, **phaselist=NULL;
#endif

partype quarklat(partype z, int *phase, latype p){
#if HINDEX==2
  const int phaselist[ORDER_GROUP][2]={{0,0}, {6,6}, {0,4}, {6,2}, 
				 {5,7}, {7,1}, {5,3}, {3,1}};
#endif
  *phase += phaselist[p][(z>>HSHIFT)&1];
  if(flip[p])z ^= 1 << HSHIFT;
  return z;
}

partype antiquarklat(partype z, int *phase, latype p){
#if HINDEX==2
  const int phaselist[ORDER_GROUP][2]={{0,0}, {2,2}, {4,0}, {6,2}, 
				 {1,3}, {7,1}, {5,3}, {7,5}};
#endif
  *phase += phaselist[p][(z>>HSHIFT)&1];
  if(flip[p])z ^= 1 << HSHIFT;
  return z;
}

/*  
    lattice group element acting on a gluon field.
    it is not very fast
 */

#define GLUE_DEBUG 0 /*  this routine is for debugging, use statelat */
#if GLUE_DEBUG
partype gluelat(partype z, latype p){
 int i,hlist[HINDEX],h2[HINDEX];
  const int *perm=permutations+(p>>HINDEX)*HINDEX;
  (void) gluedecode(z,hlist);
  for(i=0; i<HINDEX; i++)h2[i]=hlist[perm[i]];
  for(i=0; i<HINDEX; i++)if(p&(1<<i))h2[i]*=-1;
  return glueencode(z,h2);
}
#endif

/* 
    lattice group element acting on a whole state 
    *in and *out can be the same array.  

    It returns phase where the phase
    factor is equal to i^(phase/2).

    p is an element of the full lattice group.  If a subgroup
    is called, then one uses a construct of the form:

         const latlist *subg=latgroups+subgroup[multiplet];
         for(g=0; g<subg->length; g++)
 	   statelat(gc[g],ipart,type,npi,subg->e[g]);

    Since gluelat() is a BIG time user in generating bases,
    this is done "in line" Use GLUE_DEBUG to check this.
*/

int statelat(partype *out, partype *in, enum state_type type, 
	      int np, latype p){
  int i, phase=0;
#if HINDEX==1
  const partype hlist[ORDER_GROUP][2*HINDEX]={{0,1},{1,0}};
#elif HINDEX==2
  const partype hlist[ORDER_GROUP][2*HINDEX]={
    {0,1,2,3},{1,0,2,3},{0,1,3,2},{1,0,3,2},
    {2,3,0,1},{2,3,1,0},{3,2,0,1},{3,2,1,0}};
#else
  const partype **hlist=NULL;
#endif

  /* for HEAVY_LIGHT, assume that the heavy particle is a quark */
  if(type==TYPE_HEAVY_LIGHT)quarklat(0,&phase,p);

  for(i=0; i<np; i++){
    if(i==0 && type==TYPE_MESON)
      out[i]=quarklat(in[i],&phase,p);
    else if(i==np-1 && (type==TYPE_MESON || type==TYPE_HEAVY_LIGHT))
      out[i]=antiquarklat(in[i],&phase,p);
    else
#if GLUE_DEBUG  /* this is the slow way */
      out[i]=gluelat(in[i],p);
#else /* this is the fast way */
      out[i]=(in[i]&KMASK)|(hlist[p][in[i]>>HSHIFT]<<HSHIFT);
#endif      
  }

  return phase;
}

/*
   Orientation reversal of a state.
   By definition, P_1 P_2 ... orientation  = charge_conjugation.
   (The order here affects the fermion phase factor.)
   *in and *out can point to the same array.

   It returns phase where the phase factor is equal to i^(phase/2).

   for TYPE_HEAVY_LIGHT, the state returned does not
   correspond to any defined type.
*/

int orientation(partype *out, partype *in, enum state_type type, int np){
  int i, phase=0;
  partype temp;
  for(i=0; i<(np+1)/2; i++){  /* need to copy when i==np-i-1 */
    temp=in[i];           /* need temp when in[] and out[] are same */
    out[i]=in[np-i-1];
    out[np-i-1]=temp;
  }
  if(type==TYPE_MESON){
    phase += 4; /* phase from fermion anticommutator */
    /* reflection in all transverse directions */
    out[0] = quarklat(out[0],&phase,(1<<HINDEX)-1); 
    out[np-1] = antiquarklat(out[np-1],&phase,(1<<HINDEX)-1);
  } else if(type==TYPE_HEAVY_LIGHT){
    phase += 4; /* phase from fermion anticommutator */
    quarklat(0,&phase,(1<<HINDEX)-1); /* assume heavy particle is a quark */
    out[0] = quarklat(out[0],&phase,(1<<HINDEX)-1); 
  }
/* 
    Take the complex conjugate (minus the phase) because 
    orientation = P_1^* P_2^* charge_conjugation.
    *** I hope I got this part right *** 
*/
  return -phase;
}

/*
   Charge conjugation of a state.
   By definition, P_1 P_2 ... orientation  = charge_conjugation.
   (The order here affects the fermion phase factor.)
   *in and *out can point to the same array.

   It returns phase where the phase
   factor is equal to i^(phase/2).

   statelat is not called since the transverse reflection is
   especially fast to implement

   for TYPE_HEAVY_LIGHT, the state returned does not
   correspond to any defined type.
*/

int charge_conjugation(partype *out, partype *in, enum state_type type, int np){
  int phase=0,i;
  partype temp;

  for(i=0; i<(np+1)/2; i++){  /* need to copy when i==np-i-1 */
    temp=in[i];           /* need temp when in[] and out[] are same */
    out[i]=in[np-i-1];
    out[np-i-1]=temp;
  }

  /* phase from fermion anticommutator 
     assume HEAVY_LIGHT has a heavy fermion */
  if(type==TYPE_MESON||type==TYPE_HEAVY_LIGHT)phase += 4; 

  for(i=0; i<np; i++){
    if(i==0 && (type==TYPE_MESON || type==TYPE_HEAVY_LIGHT))continue;
    if(i==np-1 && type==TYPE_MESON)continue;
    out[i]^=(1<<HSHIFT);
  }

  return phase;
}

/* 
   function glueencode takes a momentum and list containing the 
   helicity and gives encoded link field.
   values of helicity are:  {1,0} {0,1} {-1,0} {0,-1}

   k2 is the discretized momentum, rounded down.

*/

partype glueencode(partype k2, int *hlist){
  int i;
  for(i=0; i<HINDEX; i++)
    if(hlist[i])break;
  return (i<<(HSHIFT+1))+((hlist[i]>0)<<HSHIFT)+(k2 & KMASK);
}

/*
  gluedecode is the inverse of glueencode.
  This routine does not shift momentum for antiperiodic boundary conditions

*/

#ifndef gluedecode /* There is also a macro version defined in basis.h */
partype gluedecode(partype x,int *hlist){
#if HINDEX==1
  hlist[0]=
#else
#if HINDEX==2
  hlist[0]=hlist[1]=0;
#else
  memset(hlist,0,HINDEX*sizeof(int));
#endif
  hlist[x>>(HSHIFT+1)] = 
#endif
    x&(1<<HSHIFT)?1:-1;
  /* return the momentum part of x */
  return x & KMASK;
}
#endif

/*
  Return a quark or given a spin {0 is spin up and 1 is spin down} 
  k2 is the quark momentum rounded down.
*/

partype quarkencode(int spin, unsigned int k2){
  return (spin<<HSHIFT) + (k2 & KMASK);
}

/*
  return quark momentum and spin given a particle encoding 
  This is the inverse of quarkencode().
*/

partype quarkdecode(partype p, int *spin){
  *spin=(p>>HSHIFT)&1;
  return p&KMASK;
}

/*
    Print out a basis state
    Here, np can be the truncation
*/

void printstate(partype *state, enum state_type type, int np){
  int j,k,spin;
  for(j=0; j<np; j++){ 
    if(state[j]==HOLE)break;
    if((type==TYPE_MESON && (j==0 || j==np-1 || state[j+1]==HOLE))
       || (type==TYPE_HEAVY_LIGHT && (j==np-1 || state[j+1]==HOLE))){
#if !QUARK_PERIODIC
      k=quarkdecode(state[j],&spin);
      if(state[j]>>HSHIFT>1)
	printf(" {%i,%2i/2}",state[j]>>HSHIFT,2*k+1);
      else
	printf(" {%c,%2i/2}",spin?'-':'+',2*k+1);
#else
      fprintf(__FILE__":  Can't print periodic quarks\n");
#endif
    } else {
      /* This expression should agree with the output of gluedecode() */
      k=state[j]&KMASK;
      spin=state[j]>>HSHIFT;
      printf(" (%2i,%i%s)",(spin&1)?+((spin>>1)+1):-((spin>>1)+1),
	   GLUE_PERIODIC?k:2*k+1,GLUE_PERIODIC?"":"/2");
    }
  }
}

/* string containing state type */

char *state_type_string(enum state_type t){
  static char string[4][20]={"glueball","meson","source","heavy-light"};
  return string[t];
}