#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "include.h"
#include "group.h"
#include "row.h"

/* 
   This file contains routines to find the remaining
   states in a basis that might have an interaction with
   a given state.

   the routines *_emerge(...) find the possible particles
   that might emerge from an interaction.  Right now, they
   are computed "on the fly"  a possible speed improvement would
   be to use look-up tables.
*/

/*
  given adjacent gluons, find possible gluons that
  may emerge from a general interaction
*/

void gggg_emerge(partype *p, int nin, interact *x){
  int i;
  partype *sp=NULL;
  partype kt=0; /* kt is the total K, rounded down */
  for(i=0; i<nin; i++)kt+=(p[i]&KMASK);
  x->np=4-nin;
  if(nin==3){
#if !GLUE_PERIODIC
    kt+=1;
#endif
    if(H_MINUS(p[0],p[1])){
      row_memory(x,1);
      x->s[0]=kt|(p[2]&(~KMASK));
    } else if(H_MINUS(p[1],p[2])){
      row_memory(x,1);
      x->s[0]=kt|(p[0]&(~KMASK));
    } else if(H_MINUS(p[0],p[2])){
      row_memory(x,1);
      x->s[0]=kt|(p[1]&(~KMASK));
    } else 
      row_memory(x,0);
  } else if(nin==2){
    partype k;
    partype h1=p[0]&(~KMASK), h2=p[1]&(~KMASK);
    if(h1==h2){
      row_memory(x,kt+1);
      for(k=0, sp=x->s; k<=kt; k++){
	sp[0]=k|h1;
	sp[1]=(kt-k)|h2;
	sp+=2;
      }
    } else if((h1^h2) == 1<<HSHIFT){
      partype i;
      row_memory(x,(kt+1)*2*HINDEX);
      for(i=0, sp=x->s; i<2*HINDEX; i++)
	for(k=0; k<=kt; k++){
	  sp[0]=k|(i<<HSHIFT);
	  sp[1]=(kt-k)|((i^1)<<HSHIFT);
	  sp+=2;
	}
    } else {
      /* orthogonal case */
      row_memory(x,(kt+1)*2);
      for(k=0, sp=x->s; k<=kt; k++){
	sp[0]=k|h1;
	sp[1]=(kt-k)|h2;
	sp[2]=k|h2;
	sp[3]=(kt-k)|h1;
	sp+=4;
      }
    }
  } else if(nin==1){
    unsigned int count=0;
    partype k1,k2,k3,h;
    partype h1=p[0]&(~KMASK);
#if !GLUE_PERIODIC
    if(kt==0){
      row_memory(x,0);
      return;
    } else
      kt-=1;
#endif
    row_memory(x,3*(2*HINDEX-1)*(kt+1)*(kt+2)/2);
    for(k1=0, sp=x->s; k1<=kt; k1++)
      for(k2=0; k2<=kt-k1; k2++){
	k3=kt-k1-k2;
	for(h=0; h<2*HINDEX; h++){
	  if(h1==(h<<HSHIFT))continue; /* Don't overcount states */
	  sp[0]=k1|h1;              
	  sp[1]=k2|((h^1)<<HSHIFT);       
	  sp[2]=k3|(h<<HSHIFT);
	  sp[3]=k1|(h<<HSHIFT);     
	  sp[4]=k2|h1;          
	  sp[5]=k3|((h^1)<<HSHIFT);
	  sp[6]=k1|((h^1)<<HSHIFT); 
	  sp[7]=k2|(h<<HSHIFT); 
	  sp[8]=k3|h1;
	  sp+=9;
	  count+=3;
	  if(count>x->length)
	    fprintf(stderr,__FILE__":  count=%i > length=%i\n",
		    count,x->length);
	}  
      }
  }
}

/*
  given adjacent quarks, find possible quarks that
  may emerge from a general interaction
*/

void qqqq_emerge(partype *p, int nin, interact *x){
  int i;
  partype *sp=NULL;
  partype kt=0;  /* total K (each rounded down) */
  partype k,h1,h2;
  if(nin!=2){
    fprintf(stderr,__FILE__":  qqq_emerge error\n");
    exit(5);
  }
  for(i=0; i<nin; i++)kt+=p[i]&KMASK;
  x->np=2;
  row_memory(x,4*(kt+1));
  for(h1=0, sp=x->s; h1<2; h1++)
    for(h2=0; h2<2; h2++)
      for(k=0; k<=kt; k++){
	sp[0]=k|(h1<<HSHIFT);
	sp[1]=(kt-k)|(h2<<HSHIFT);
	sp+=2;
      }
}

/*
  given adjacent quark and source, find possible quarks that
  may emerge from a general interaction

  kt_left is twice the K cut-off
*/

void ssqq_emerge(partype *p, int nin, int kt_left, interact *x){
  partype *sp=NULL;
  partype kt=0; /* K, rounded down */
  partype k,h1;
  if(nin!=1){
    fprintf(stderr,__FILE__":  ssqq_emerge error\n");
    exit(5);
  }
  x->np=1;
  row_memory(x,kt+1);
  for(h1=0, sp=x->s; h1<2; h1++)
    for(k=0; 2*k+1<=kt_left; k++){
      sp[0]=(h1<<HSHIFT)|k;
      sp++;
    }
}

/*
       for a general qqgg interaction, find the particles
       that can emerge. 
       p[0] is the quark
*/

void qqgg_emerge(partype *p, int nin, interact *x){
  partype k,i,j;
  partype *sp=NULL;
  partype kt=0;  /* kt is total K, each rounded down */
  for(i=0; i<nin; i++)kt+=p[i]&KMASK; 
  x->np=4-nin;
  if(nin==3){
#if GLUE_PERIODIC
    row_memory(x,2);
    x->s[0]=kt;
    x->s[1]=kt|(1<<HSHIFT);
#endif
  }else if(nin==2){
    row_memory(x,4*HINDEX*(kt+1));
    for(i=0, sp=x->s; i<2; i++)
      for(j=0; j<2*HINDEX; j++)
	for(k=0; k<=kt; k++){
	  sp[0]=k|(i<<HSHIFT);
	  sp[1]=(kt-k)|(j<<HSHIFT);
	  sp+=2;
	}
  } else if(nin==1) {
    unsigned int count=0;
    partype k1,k2,k3,h1,h2,j;
    row_memory(x,2*(2*HINDEX)*(2*HINDEX)*(kt+1)*(kt+2)/2);
    for(k1=0, sp=x->s; k1<=kt; k1++)
      for(k2=0; k2<=kt-k1; k2++){
	k3=kt-k1-k2;
	for(j=0; j<2; j++)
	  for(h1=0; h1<2*HINDEX; h1++)
	    for(h2=0; h2<2*HINDEX; h2++){
	      sp[0]=k1|(j<<HSHIFT);    
	      sp[1]=k2|(h1<<HSHIFT); 
	      sp[2]=k3|(h2<<HSHIFT); 
	      sp+=3;
	      count++;
	      if(count>x->length)
		fprintf(stderr,__FILE__":  count=%i > length=%i\n",
			count,x->length);
	    }
      }
  }
}

/*
      ssgg is gives a list of particles that can
      emerge from a static source-gluon interation.

      kt_left twice the K cutoff minus the total light field momenta.

      npin is the number of incoming link fields
*/

void ssgg_emerge(partype *p, int nin, int kt_left, interact *x){
  int i;
  partype *sp=NULL;
  for(i=0; i<nin; i++)kt_left+=2*(p[i]&KMASK)+(GLUE_PERIODIC?0:1);
  x->np=2-nin;
  if(nin==2){
    if(((p[0]^p[1])&(~KMASK))==1<<HSHIFT)
      row_memory(x,1);
    else
      row_memory(x,0);
  } else if(nin==1){
    partype k;
    partype h1=p[0]&(~KMASK);
    row_memory(x,(kt_left+(GLUE_PERIODIC?2:1))/2);
    for(k=0, sp=x->s; 2*k+(GLUE_PERIODIC?0:1)<=kt_left; k++){
      sp[0]=h1|k;
      sp++;
    }
  } else if(nin==0){
    partype k1,k2,h;
    row_memory(x,2*HINDEX*(kt_left/2+(GLUE_PERIODIC?1:0))*
      (kt_left/2+1+(GLUE_PERIODIC?1:0))/2);
    for(k1=0, sp=x->s; k1<=kt_left/2; k1++)
      for(k2=0; k1+k2+(GLUE_PERIODIC?0:1)<=kt_left/2; k2++){
	for(h=0; h<2*HINDEX; h++){
	  sp[0]=(h<<HSHIFT)|k1;    
	  sp[1]=((h^1)<<HSHIFT)|k2;  
	  sp+=2;
	}
      }
  }
}

/*
       for a general qqg interaction, find the particles
       that can emerge. 
       p[0] is the quark
*/

void qqg_emerge(partype *p, int nin, interact *x){
  int i;
  partype *sp=NULL;
  partype kt=0;   /* kt is the total K, each rounded down */
  for(i=0; i<nin; i++)kt+=p[i]&KMASK;
  x->np=3-nin;
  if(nin==2){
#if GLUE_PERIODIC
    row_memory(x,2);
    x->s[0]=kt;
    x->s[1]=kt|(1<<HSHIFT);
#endif
  } else if(nin==1){
    partype k1,j,h;
    row_memory(x,2*2*HINDEX*(kt+1));
    for(k1=0, sp=x->s; k1<=kt; k1++)
      for(j=0; j<2; j++)
	for(h=0; h<2*HINDEX; h++){
	  sp[0]=k1|(j<<HSHIFT);    
	  sp[1]=(kt-k1)|(h<<HSHIFT); 
	  sp+=2;
	}    
  }
}
 
/**************************************************************/

/* 
    row_*(...) finds the remaining
    states in the basis bb having a possible nonzero matrix
    element with state x and puts them in y.
    It returns states in the basis occuring after or equal to x.
    It returns the resulting list of locations.  
   
    The calling routine must allocate memory for y.

    Throughout, the following relation should hold:

    if(npout!=npx-nin+s.np) 
      fprintf(stderr,__FILE__":  interaction disaster\n");

*/

#define FOUR 4  /* number of particles in link interactions */
#define THREE 3 /* number of particles in Yukawa 3-point interactions */
#define TWO 2   /* number of particles in source interactions */

/*
  Fill remaining slots with holes, well-order state,
  check if state is valid,
  and add state to sorted basis if it occurs *after* 
  the first state in the basis.
*/

void fix_add(cbase *y, full_basis *bb, partype *z, int np){
  int l;
#if 0       /* debugging print */
  printf("    final:");printstate(z,bb->type,np);printf("\n");
#endif
  if(!zeromodetest(z,bb->type,np))return; /* check for too many zero modes */
  /* check for non-zero norm 
     This is a large time user for constructing a
     matrix.  Perhaps construct a list of zero norm states
     and perform a search. Or check for zero norm only if state is
     not found in basis. ?????? */
  if(0==innerproduct(z,z,bb->type,np,bb->charge,bb->o,bb->multi))return;
  wellorder(z,bb->type,np,bb->charge,bb->o,subgroup[bb->multi]);
  for(l=np; l<bb->nptrunc; l++)z[l]=HOLE;
  for(l=0; l<bb->nptrunc; l++){
    if(z[l]<y->s[l])return;
    else if(z[l]>y->s[l]){
      y->length=addstate(y->s,z,y->length,bb->nptrunc);
      return;
    }
  }
}

void row_glue(cbase *y, full_basis *bb, partype *x, int npx){
  int nin,npout;
  int i,j,l;
  partype z[NPMAX];
  partype xx[2*NPMAX];
  interact s=INTERACT_INIT; /* object holding results of interactions */ 
  y->length=0;
  y->length=addstate(y->s,x,y->length,bb->nptrunc);
  if(bb->type!=TYPE_GLUE){
    fprintf(stderr,__FILE__":  Wrong type of interaction");
    exit(55);
  }
  memcpy(xx,x,npx*sizeof(partype));
  memcpy(xx+npx,x,npx*sizeof(partype));
  for(nin=1; nin<FOUR; nin++){
    npout=npx+FOUR-2*nin;
    if(nin>npx || npout>bb->nptrunc)continue;
    for(i=0; i<npx; i++){
#if 0        /* debugging print */
      printf("Starting  ");
      printstate(xx+i,bb->type,npx);
      printf(" nin=%i\n",nin);
#endif
      gggg_emerge(xx+i,nin,&s);
      for(j=0; j<s.length; j++){
	for(l=0; l<s.np; l++)
	  z[l]=s.s[s.np*j+l];
	for(l=0; l<npx-nin; l++)
	  z[s.np+l]=xx[l+i+nin];
	fix_add(y,bb,z,npout);
      }
    }
  }

  free_row_memory(&s);
}



void row_source(cbase *y, full_basis *bb, partype *x, int npx){
  int nin,npout;
  int i,j,l;
  partype z[NPMAX];
  partype kt=bb->kt; /* twice the total K cutoff minus light field total K */
  partype xx[NPMAX];
  interact s=INTERACT_INIT; /* object holding results of interactions */ 
  y->length=0;
  y->length=addstate(y->s,x,y->length,bb->nptrunc);
  if(bb->type!=TYPE_SOURCE){
    fprintf(stderr,__FILE__":  Wrong type of interaction");
    exit(55);
  }
  for(i=0; i<npx; i++)xx[i]=x[npx-1-i];
  for(i=0; i<npx; i++)kt-=2*(x[i]&KMASK)+(GLUE_PERIODIC?0:1);

  /* gggg interactions acting on interior of string */

  for(nin=1; nin<FOUR; nin++){
    npout=npx+FOUR-2*nin;
    if(npout>bb->nptrunc)continue;
    for(i=0; i<=npx-nin; i++){
      gggg_emerge(x+i,nin,&s);
      for(j=0; j<s.length; j++){
	for(l=0; l<i; l++)
	  z[l]=x[l];
	for(l=0; l<s.np; l++)
	  z[l+i]=s.s[s.np*j+l];
	for(l=i; l<npx-nin; l++)
	  z[s.np+l]=x[l+nin];
	fix_add(y,bb,z,npout);
      }
    }
  }

  /* ssgg interactions */  

  for(nin=0; nin<=TWO; nin++){
    npout=npx+TWO-2*nin;
    if(nin>npx || npout>bb->nptrunc)continue;
#if 0        /* debugging print */
    printf("Starting  ");
    printstate(x,bb->type,npx);
    printf(" nin=%i kt=%i\n",nin,kt);
#endif
    ssgg_emerge(x,nin,kt,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=x[l+nin];
      fix_add(y,bb,z,npout);
    }
    ssgg_emerge(xx,nin,kt,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<npx-nin; l++)
	z[l]=x[l];
      for(l=0; l<s.np; l++)
	z[npout-1-l]=s.s[s.np*j+l];
      fix_add(y,bb,z,npout);
    }
  }

  free_row_memory(&s);
}

void row_meson(cbase *y, full_basis *bb, partype *x, int npx){
  int nin,npout;
  int i,j,l;
  partype z[NPMAX];
  partype xx[NPMAX];
  interact s=INTERACT_INIT; /* object holding results of interactions */ 
  y->length=0;
  y->length=addstate(y->s,x,y->length,bb->nptrunc);
  if(bb->type!=TYPE_MESON){
    fprintf(stderr,__FILE__":  Wrong type of interaction");
    exit(55);
  }

  /* qqqq interactions */

  if(npx==2){ 
#if 0        /* debugging print */
    printf("Starting qqqq:  ");
    printstate(x,bb->type,npx);
    printf("\n");
#endif
    qqqq_emerge(x,npx,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      fix_add(y,bb,z,s.np);
    }
  }

  /* gggg gauge interactions */

  for(nin=1; nin<FOUR; nin++){
    npout=npx+FOUR-2*nin;
    if(nin>npx-2 || npout>bb->nptrunc)continue;
    for(i=1; i<npx-nin; i++){
#if 0        /* debugging print */
      printf("Starting gggg:");
      printstate(x+i,bb->type,npx);
      printf(" nin=%i\n",nin);
#endif
      gggg_emerge(x+i,nin,&s);
      for(j=0; j<s.length; j++){
	for(l=0; l<i; l++)
	  z[l]=x[l];
	for(l=0; l<s.np; l++)
	  z[l+i]=s.s[s.np*j+l];
	for(l=i; l<npx-nin; l++)
	  z[s.np+l]=x[l+nin];
	fix_add(y,bb,z,npout);
      }
    }
  }
  
  charge_conjugation(xx,x,bb->type,npx);

  /* qqg Yukawa interactions */
  
  for(nin=1; nin<THREE; nin++){
    npout=npx+THREE-2*nin;
    if(nin>npx-1 || npout>bb->nptrunc)continue;
#if 0        /* debugging print */
    printf("Starting qqg:  ");
    printstate(x,bb->type,npx);
    printf(" nin=%i\n",nin);
#endif
   qqg_emerge(x,nin,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=x[l+nin];
      fix_add(y,bb,z,npout);
    }
    qqg_emerge(xx,nin,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=xx[l+nin];
      charge_conjugation(z,z,bb->type,npout);
      fix_add(y,bb,z,npout);
    }
  }

  /*  qqgg Yukawa and gauge interactions    */  

  for(nin=1; nin<FOUR; nin++){
    npout=npx+FOUR-2*nin;
#if 0        /* debugging print */
    printf("Starting qqgg:  ");
    printstate(x,bb->type,npx);
    printf(" nin=%i\n",nin);
#endif
    if(nin>npx-1 || npout>bb->nptrunc)continue;
    qqgg_emerge(x,nin,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=x[l+nin];
      fix_add(y,bb,z,npout);
    }
    qqgg_emerge(xx,nin,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=xx[l+nin];
      charge_conjugation(z,z,bb->type,npout);
      fix_add(y,bb,z,npout);
    }
  }

  free_row_memory(&s);
}

void row_heavy_light(cbase *y, full_basis *bb, partype *x, int npx){
  int nin,npout;
  int i,j,l;
  partype z[NPMAX];
  partype kt=bb->kt; /* twice K cutoff minus light field total K */
  partype xx[NPMAX];
  interact s=INTERACT_INIT; /* object holding results of interactions */ 
  y->length=0;
  y->length=addstate(y->s,x,y->length,bb->nptrunc);
  if(bb->type!=TYPE_HEAVY_LIGHT){
    fprintf(stderr,__FILE__":  Wrong type of interaction");
    exit(55);
  }

  for(i=0; i<npx-1; i++)kt-=2*(x[i]&KMASK)+(GLUE_PERIODIC?0:1);
  kt-=2*(x[npx-1]&KMASK)+(QUARK_PERIODIC?0:1);

  /* ssqq interactions */

  if(npx==1){ 
#if 0        /* debugging print */
    printf("Starting qqqq:  ");
    printstate(x,bb->type,npx);
    printf("\n");
#endif
    ssqq_emerge(x,npx,bb->kt,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      fix_add(y,bb,z,s.np);
    }
  }

  /* gggg gauge interactions */

  for(nin=0; nin<FOUR; nin++){
    npout=npx+FOUR-2*nin;
    if(npout>bb->nptrunc)continue;
    for(i=0; i<npx-nin; i++){
#if 0        /* debugging print */
      printf("Starting gggg:");
      printstate(x+i,bb->type,npx);
      printf(" nin=%i\n",nin);
#endif
      gggg_emerge(x+i,nin,&s);
      for(j=0; j<s.length; j++){
	for(l=0; l<i; l++)
	  z[l]=x[l];
	for(l=0; l<s.np; l++)
	  z[l+i]=s.s[s.np*j+l];
	for(l=i; l<npx-nin; l++)
	  z[s.np+l]=x[l+nin];
	fix_add(y,bb,z,npout);
      }
    }
  }
  
  charge_conjugation(xx,x,bb->type,npx);

  /* qqg Yukawa interactions */
  
  for(nin=0; nin<THREE; nin++){
    npout=npx+THREE-2*nin;
    if(nin>npx || npout>bb->nptrunc)continue;
#if 0        /* debugging print */
    printf("Starting qqg:  ");
    printstate(x,bb->type,npx);
    printf(" nin=%i\n",nin);
#endif
    qqg_emerge(xx,nin,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=xx[l+nin];
      charge_conjugation(z,z,bb->type,npout);
      fix_add(y,bb,z,npout);
    }
  }

  /*  qqgg Yukawa and gauge interactions    */  

  for(nin=0; nin<FOUR; nin++){
    npout=npx+FOUR-2*nin;
#if 0        /* debugging print */
    printf("Starting qqgg:  ");
    printstate(x,bb->type,npx);
    printf(" nin=%i\n",nin);
#endif
    if(nin>npx || npout>bb->nptrunc)continue;
    qqgg_emerge(xx,nin,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=xx[l+nin];
      charge_conjugation(z,z,bb->type,npout);
      fix_add(y,bb,z,npout);
    }
  }

  for(nin=0; nin<=TWO; nin++){
    npout=npx+TWO-2*nin;
    if(nin>npx || npout>bb->nptrunc)continue;
    ssgg_emerge(x,nin,kt,&s);
    for(j=0; j<s.length; j++){
      for(l=0; l<s.np; l++)
	z[l]=s.s[s.np*j+l];
      for(l=0; l<npx-nin; l++)
	z[s.np+l]=x[l+nin];
	fix_add(y,bb,z,npout);
      }
    }

  free_row_memory(&s);
}


void row_memory(interact *y, size_t length){
  size_t temp_length=length*y->np*sizeof(partype);
  if(temp_length>y->s_length){
    y->s=realloc(y->s,y->s_length=temp_length);
    if(y->s==NULL){
      fprintf(stderr,__FILE__":  Can't allocate memory in row_memory\n");
      exit(11);
    }
  }
  y->length=length;
}

void free_row_memory(interact *s){
  free(s->s);s->s=NULL;
  s->s_length=0;
}