#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "improved.h"
#include "interact.h"

/* 
  There are two ways that the code has been speeded up:
  first, use of the mod operator "%" has been removed.
  second, "memcmp" is slower than in-line code
*/

/*
  Initialize various parameters for interaction calculations.
  This routine should eventually take over the various other
  initialization routines.
*/

void init_params(int nptrunc, unsigned int kt, int momflag, element *apperp){
  int i;
  working.nptrunc=nptrunc;
  working.kt=kt;
  working.momflag=momflag;
  for(i=0; i<HINDEX; i++)working.apperp[i]=apperp[i];
}

/*
  Below are subroutines for glueball states.  They should
  work for 2+1 and 3+1 dimensions.

  multi is a basis of irreducable representations of the
  lattice symmetries 

  Define lists of phases used by various subroutines
  Half integers are defined for
  phaselist and phaselistv when HINDEX>1.
*/

static complex_element phaselist[2*KKMAX][HINDEX];
static element phaseangleu[HINDEX];
static complex_element phaselistv[2*NPMAX][HINDEX];

void makephaselist(element *angle,int kt){
  int i,j;
  element theta;
  for(j=0; j<HINDEX; j++){
    theta=0.5*2.0*angle[j]/(element) kt;
    for(i=0; i<2*KKMAX; i++){
      phaselist[i][j].re=cos(theta*(element) i);
      phaselist[i][j].im=sin(theta*(element) i);
    }
    phaseangleu[j]=angle[j];
    for(i=0; i<2*NPMAX; i++){
      phaselistv[i][j].re=cos(0.5*angle[j]*(element) i);
      phaselistv[i][j].im=sin(0.5*angle[j]*(element) i);
    }
  }
#if 0  /* debug print */
  printf("phases %f %f %f %f\n",phaselist[0].re,phaselist[0].im,
		phaselist[1].re,phaselist[1].im); 
#endif
}

/*
  find phase factor for given transition, all link fields.

  This could be further speeded up by using a look-up table
  to do the sines and cosines.
*/

void glue_phase(complex_element *z, partype *left, int npl,
		partype *right, int npr){
  int i,k;
  element x[HINDEX],phase=0;
  for(k=0; k<HINDEX; k++)x[k]=0;
  for(i=0; i<npl; i++){
    for(k=0; k<i; k++)
      x[left[k]>>(HSHIFT+1)]+=GE(left[i])*(left[k]&(1<<HSHIFT)?-1:1);
    x[left[i]>>(HSHIFT+1)]+=GE(left[i])*(left[i]&(1<<HSHIFT)?-0.5:0.5);
  }
  for(i=0; i<npr; i++){
    for(k=0; k<i; k++)
      x[right[k]>>(HSHIFT+1)]-=GE(right[i])*(right[k]&(1<<HSHIFT)?-1:1);
    x[right[i]>>(HSHIFT+1)]-=GE(right[i])*(right[i]&(1<<HSHIFT)?-0.5:0.5);
  }
  for(i=0; i<HINDEX; i++)phase+=working.apperp[i]*x[i];
  phase/=0.5*(element) working.kt; /* working.kt is twice the momentum */
  z->re=cos(phase);
  z->im=sin(phase);
}

/*
    Find matrix element of the Hamiltonian assuming cyclic permutaitons
    only for the states. 

    It might help to "unroll" the innermost loops since
    they have fixed iteraction number.  However, no real
    improvement was seen with gcc compiler.

    This is the interaction for the case of zero transverse
    momentum.
*/

void interact_glueball(element *x, partype *left, int npl, partype *right, 
		 int npr, element *couplings){
  partype l1,l2,l3,r1,r2,r3;
  partype ka,kb;
  fun22_forms *row;
  fun13_forms *sow;
  int hl1[HINDEX],hl2[HINDEX],hl3[HINDEX],hr1[HINDEX],hr2[HINDEX],hr3[HINDEX];
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
  int il,ir;
  int j;
  memcpy(leftleft,left,npl*sizeof(partype));
  memcpy(leftleft+npl,left,npl*sizeof(partype));
  memcpy(rightright,right,npr*sizeof(partype));
  memcpy(rightright+npr,right,npr*sizeof(partype));
  *x=0.0;
  if(npl==npr && npl >1){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
	  l1=gluedecode(leftleft[il],hl1); l2=gluedecode(leftleft[il+1],hl2);
	  r1=gluedecode(rightright[ir],hr1); r2=gluedecode(rightright[ir+1],hr2);
	  sort_k(&ka,&kb,l1,l2,r1,r2);
	  row=&fun22[l1+l2][ka][kb];
          if(EQUAL(hl1,hr1)){
	    *x+=couplings[MM]*row->mass; 
            *x+=couplings[GGN]*row->ggn;
            if(EQUAL(hl1,hl2))
	      *x+=couplings[LAMBDA_2]*row->lambda;
 /* pinching terms are equal to the half and 1 the lambda_1 and lambda_4 
    interactions, respectively, for npl=npr=2*/
	    else if(MINUS(hl1,hl2)){
	      *x+=couplings[LAMBDA_1]*row->lambda;
	      if(npl==2)*x+=0.5*couplings[LAMBDA_3]*row->lambda;
	    } else {
	      *x+=couplings[LAMBDA_4]*row->lambda;
	      if(npl==2 && MINUS(hl1,hl2)) /* no lambda_5 for winding modes. */
		*x+=couplings[LAMBDA_5]*row->lambda;
	    } 
	  }
          if(MINUS(hl1,hl2)){
            if(npl!=2)
	      *x-=couplings[GGN]*row->annihilate;
            if(EQUAL(hl1,hr1)){
	      *x+=couplings[LAMBDA_1]*row->lambda;
	      if(npl==2)*x+=0.5*couplings[LAMBDA_3]*row->lambda;
	    } else if(MINUS(hl1,hr1))
	      *x+=couplings[LAMBDA_2]*row->lambda;
	    else {
	      *x+=couplings[LAMBDA_4]*row->lambda;
	      if(npl==2 && MINUS(hl1,hl2)) /* lambda_5 not for winding modes. */
		*x+=couplings[LAMBDA_5]*row->lambda;
	    }
          } 
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hl2))
	    *x-=couplings[BETA]*row->lambda;
        } 
      } 
  } else if(npl+2==npr){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
	  l1=gluedecode(leftleft[il],hl1); r1=gluedecode(rightright[ir],hr1); 
	  r2=gluedecode(rightright[ir+1],hr2); r3=gluedecode(rightright[ir+2],hr3);
          if(EQUAL(hl1,hr3)){
	    sow=&fun13[l1][r1][r2];
	    *x+=couplings[GGN]*sow->ggn;
            if(EQUAL(hl1,hr2))
	      *x+=couplings[LAMBDA_2]*sow->lambda;
	    else if(MINUS(hl1,hr2))
	      *x+=couplings[LAMBDA_1]*sow->lambda;
	    else
	      *x+=couplings[LAMBDA_4]*sow->lambda;
	  }
          if(EQUAL(hl1,hr1)){
	    sow=&fun13[l1][r3][r2];
	    *x+=couplings[GGN]*sow->ggn;
            if(EQUAL(hl1,hr2))
	      *x+=couplings[LAMBDA_2]*sow->lambda;
	    else if(MINUS(hl1,hr2))
	      *x+=couplings[LAMBDA_1]*sow->lambda;
	    else
	      *x+=couplings[LAMBDA_4]*sow->lambda;
	  }
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hr1)){
	    sow=&fun13[l1][r1][r2];
	    *x-=couplings[BETA]*sow->lambda;
	  }
	} 
      } 
  }
  else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
	  l1=gluedecode(leftleft[il],hl1); l2=gluedecode(leftleft[il+1],hl2);
	  l3=gluedecode(leftleft[il+2],hl3);r1=gluedecode(rightright[ir],hr1);
          if(EQUAL(hr1,hl3)){
	    sow=&fun13[r1][l1][l2];
	    *x+=couplings[GGN]*sow->ggn;
            if(EQUAL(hr1,hl2))
	      *x+=couplings[LAMBDA_2]*sow->lambda;
	    else if(MINUS(hr1,hl2))
	      *x+=couplings[LAMBDA_1]*sow->lambda;
	    else
	      *x+=couplings[LAMBDA_4]*sow->lambda;
	  }
          if(EQUAL(hr1,hl1)){
	    sow=&fun13[r1][l3][l2];
	    *x+=couplings[GGN]*sow->ggn;
            if(EQUAL(hr1,hl2))
	      *x+=couplings[LAMBDA_2]*sow->lambda;
	    else if(MINUS(hr1,hl2))
	      *x+=couplings[LAMBDA_1]*sow->lambda;
	    else
	      *x+=couplings[LAMBDA_4]*sow->lambda;
	  }
	  if(EQUAL(hr1,hl2)&&ORTHO(hl1,hr1)){
	    sow=&fun13[r1][l1][l2];
	    *x-=couplings[BETA]*sow->lambda;
	  }
	} 
      } 
  }
/* here we take the case that the matrix element is the mass */
  else if(npl==1 && npr==1){
    *x=0.5*couplings[MM]/
      ((element) gluedecode(left[0],hl1)+(GLUE_PERIODIC?0.0:0.5));
  }
  return;
}

/*
    Find matrix element of the Hamiltonian for open strings. 
    This is the interaction for the case of zero transverse
    momentum.  Unlike the glueball case, we cannot assume that 
    the winding number of the left and right states match:  this
    must be checked explicitly.

    There are no "pinching" terms and no 1-gluon mass in this case

    Here, the sum *x is not initialized.
*/

void interact_glue(element *x, partype *left, int npl, partype *right, 
		 int npr, element *couplings){
  int i,j,jj;
  partype ka,kb;
  fun22_forms *row;
  fun13_forms *sow;
  int l[NPMAX],r[NPMAX];

  for(i=0; i<npl; i++)l[i]=left[i]&KMASK;
  for(i=0; i<npr; i++)r[i]=right[i]&KMASK;

  if(npl==npr && npl>1){
    for(j=0; j<npl-1; j++){
      /* If there are any mismatches on the left 
	 hand side, then there is no hope for further matches */
      if(j>0 && left[j-1]!=right[j-1])break;
      for(jj=j+2; jj<npl; jj++)  /* 4-points */
	if(left[jj]!=right[jj])break;
      if(jj<npl)continue;  /* try next j if spectators don't match */
      sort_k(&ka,&kb,l[j],l[j+1],r[j],r[j+1]);
      row=&fun22[l[j]+l[j+1]][ka][kb];
      if(H_EQUAL(left[j],right[j]) && H_EQUAL(left[j+1],right[j+1])){
       	*x+=couplings[MM]*row->mass+couplings[GGN]*row->ggn;
	if(H_EQUAL(left[j],left[j+1]))
	  *x+=couplings[LAMBDA_2]*row->lambda;
	else if(H_MINUS(left[j],left[j+1]))
	  *x+=couplings[LAMBDA_1]*row->lambda;
	else
	  *x+=couplings[LAMBDA_4]*row->lambda;
      } 
      if(H_MINUS(left[j],left[j+1]) && H_MINUS(right[j],right[j+1])){
	/* Note that the instantaneous-annihilation graph
	   can also act in the two particle sector, unlike
	   in the closed loop case. */
	*x-=couplings[GGN]*row->annihilate;
	if(H_EQUAL(left[j],right[j]))
	  *x+=couplings[LAMBDA_1]*row->lambda;
	else if(H_MINUS(left[j],right[j]))
	  *x+=couplings[LAMBDA_2]*row->lambda;
	else
	  *x+=couplings[LAMBDA_4]*row->lambda;
      } 
      if(H_EQUAL(left[j],right[j+1]) && H_EQUAL(left[j+1],right[j]) &&
	 H_ORTHO(left[j],left[j+1]))
	*x-=couplings[BETA]*row->lambda;
    }
    /* The "pinching" terms are suppressed in the open string although
       one could introduce terms like (\phi^\d M \phi)^2  */
  }

  else if(npl>0 && npl+2==npr){
    for(j=0; j<npl; j++){
      if(j>0 && left[j-1]!=right[j-1])break; /* larger j won't work either */
      for(jj=j+1; jj<npl; jj++)
	if(left[jj]!=right[jj+2])break;
      if(jj<npl)continue;  /* go to next j if there is a mismatch */
      if(H_EQUAL(left[j],right[j+2]) && H_MINUS(right[j],right[j+1])){
	sow=&fun13[l[j]][r[j]][r[j+1]];
	*x+=couplings[GGN]*sow->ggn;
	if(H_EQUAL(left[j],right[j+1]))
	  *x+=couplings[LAMBDA_2]*sow->lambda;
	else if(H_MINUS(left[j],right[j+1]))
	  *x+=couplings[LAMBDA_1]*sow->lambda;
	else
	  *x+=couplings[LAMBDA_4]*sow->lambda;
      }
      if(H_EQUAL(left[j],right[j]) && H_MINUS(right[j+1],right[j+2])){
	sow=&fun13[l[j]][r[j+2]][r[j+1]];
	*x+=couplings[GGN]*sow->ggn;
	if(H_EQUAL(left[j],right[j+1]))
	  *x+=couplings[LAMBDA_2]*sow->lambda;
	else if(H_MINUS(left[j],right[j+1]))
	  *x+=couplings[LAMBDA_1]*sow->lambda;
	else
	  *x+=couplings[LAMBDA_4]*sow->lambda;
      }
      if(H_EQUAL(left[j],right[j+1]) && H_MINUS(right[j],right[j+2])
	 && H_ORTHO(left[j],right[j])){
	sow=&fun13[l[j]][r[j]][r[j+1]];
	*x-=couplings[BETA]*sow->lambda;
      }
    }  
  }

  else if(npr>0 && npr+2 == npl){
    for(j=0; j<npr; j++){
      if(j>0 && left[j-1]!=right[j-1])break; /* larger j won't work either */
      for(jj=j+1; jj<npr; jj++)
	if(left[jj+2]!=right[jj])break;
      if(jj<npr)continue; /* go to next j if there is a mismatch */
      if(H_EQUAL(right[j],left[j+2]) && H_MINUS(left[j],left[j+1])){
	sow=&fun13[r[j]][l[j]][l[j+1]];
	*x+=couplings[GGN]*sow->ggn;
	if(H_EQUAL(right[j],left[j+1]))
	  *x+=couplings[LAMBDA_2]*sow->lambda;
	else if(H_MINUS(right[j],left[j+1]))
	  *x+=couplings[LAMBDA_1]*sow->lambda;
	else
	  *x+=couplings[LAMBDA_4]*sow->lambda;
      }
      if(H_EQUAL(right[j],left[j]) && H_MINUS(left[j+1],left[j+2])){
	sow=&fun13[r[j]][l[j+2]][l[j+1]];
	*x+=couplings[GGN]*sow->ggn;
	if(H_EQUAL(right[j],left[j+1]))
	  *x+=couplings[LAMBDA_2]*sow->lambda;
	else if(H_MINUS(right[j],left[j+1]))
	  *x+=couplings[LAMBDA_1]*sow->lambda;
	else
	  *x+=couplings[LAMBDA_4]*sow->lambda;
      }
      if(H_EQUAL(right[j],left[j+1]) && H_MINUS(left[j],left[j+2])
	 && H_ORTHO(left[j],right[j])){
	sow=&fun13[r[j]][l[j]][l[j+1]];
	*x-=couplings[BETA]*sow->lambda;
      } 
    } 
  }

#if 1  /* debug flag */
  if(*x==0.0 || fabs(*x)<1.0e10) return;
  printf("   last x=%f:  ",*x);
  printstate(left,TYPE_SOURCE,npl);
  printf(", ");
  printstate(right,TYPE_SOURCE,npr);
  printf("\n");
#endif

  return;
}

/*
    Find matrix element of the Hamiltonian for open strings,
    nonzero transverse momentum.
    There are no "pinching" terms and no 1-gluon mass in this case

    Here, the sum *x is not initialized.
*/

void interact_gluec(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings){
  partype ka,kb;
  fun22_forms *row;
  fun13_forms *sow;
  complex_element ztemp;
  int l[NPMAX],r[NPMAX];
  int i,j,jj;

  for(i=0; i<npl; i++)l[i]=left[i]&KMASK;
  for(i=0; i<npr; i++)r[i]=right[i]&KMASK;

  if(npl==npr && npl>1){
    for(j=0; j<npl-1; j++){
      /* If there are any mismatches on the left 
	 hand side, then there is no hope for further matches */
      if(j>0 && left[j-1]!=right[j-1])break;
      for(jj=j+2; jj<npl; jj++) /* 4-points */
	if(left[jj]!=right[jj])break;
      if(jj<npl)continue;  /* try next j if spectators don't match */
      sort_k(&ka,&kb,l[j],l[j+1],r[j],r[j+1]);
      row=&fun22[l[j]+l[j+1]][ka][kb];
      glue_phase(&ztemp,left+j,2,right+j,2);
      if(H_EQUAL(left[j],right[j]) && H_EQUAL(left[j+1],right[j+1])){
	ZTIMES(x,&ztemp,couplings[MM]*row->mass+couplings[GGN]*row->ggn);
	if(H_EQUAL(left[j],left[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_2]*row->lambda);
	else if(H_MINUS(left[j],left[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_1]*row->lambda);
	else
	  ZTIMES(x,&ztemp,couplings[LAMBDA_4]*row->lambda);
      } 
      if(H_MINUS(left[j],left[j+1]) && H_MINUS(right[j],right[j+1])){
	/* Note that the instantaneous-annihilation graph
	   can also act in the two particle sector, unlike
	   in the closed loop case. */
	ZTIMES(x,&ztemp,-couplings[GGN]*row->annihilate);
	if(H_EQUAL(left[j],right[j]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_1]*row->lambda);
	else if(H_MINUS(left[j],right[j]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_2]*row->lambda);
	else
	  ZTIMES(x,&ztemp,couplings[LAMBDA_4]*row->lambda);
      } 
      if(H_EQUAL(left[j],right[j+1]) && H_EQUAL(left[j+1],right[j]) && 
		 H_ORTHO(left[j],left[j+1]))
	ZTIMES(x,&ztemp,-couplings[BETA]*row->lambda);
    }
    /* The "pinching" terms are suppressed in the open string although
       one could introduce terms like (\phi^\d M \phi)^2  */
  }

  else if(npl+2 == npr && npl>0){
    for(j=0; j<npl; j++){
      if(j>0 && left[j-1]!=right[j-1])break; /* larger j won't work either */
      for(jj=j+1; jj<npl; jj++)
	if(left[jj]!=right[jj+2])break;
      if(jj<npl)continue;  /* go to next j if there is a mismatch */
      glue_phase(&ztemp,left+j,1,right+j,3);
      if(H_EQUAL(left[j],right[j+2]) && H_MINUS(right[j],right[j+1])){
	sow=&fun13[l[j]][r[j]][r[j+1]];
	ZTIMES(x,&ztemp,couplings[GGN]*sow->ggn);
	if(H_EQUAL(left[j],right[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_2]*sow->lambda);
	else if(H_MINUS(left[j],right[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_1]*sow->lambda);
	else
	  ZTIMES(x,&ztemp,couplings[LAMBDA_4]*sow->lambda);
      }
      if(H_EQUAL(left[j],right[j]) && H_MINUS(right[j+1],right[j+2])){
	sow=&fun13[l[j]][r[j+2]][r[j+1]];
	ZTIMES(x,&ztemp,couplings[GGN]*sow->ggn);
	if(H_EQUAL(left[j],right[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_2]*sow->lambda);
	else if(H_MINUS(left[j],right[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_1]*sow->lambda);
	else
	  ZTIMES(x,&ztemp,couplings[LAMBDA_4]*sow->lambda);
      }
      if(H_EQUAL(left[j],right[j+1]) && H_MINUS(right[j],right[j+2]) && 
	 H_ORTHO(left[j],right[j])){
	sow=&fun13[l[j]][r[j]][r[j+1]];
	ZTIMES(x,&ztemp,-couplings[BETA]*sow->lambda);
      }
    }  
  }
  else if(npr+2 == npl && npr>0){
    for(j=0; j<npr; j++){
      if(j>0 && left[j-1]!=right[j-1])break; /* larger j won't work either */
      for(jj=j+1; jj<npr; jj++)
	if(left[jj+2]!=right[jj])break;
      if(jj<npr)continue; /* go to next j if there is a mismatch */
      glue_phase(&ztemp,left+j,3,right+j,1);
      if(H_EQUAL(right[j],left[j+2]) && H_MINUS(left[j],left[j+1])){
	sow=&fun13[r[j]][l[j]][l[j+1]];
	ZTIMES(x,&ztemp,couplings[GGN]*sow->ggn);
	if(H_EQUAL(right[j],left[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_2]*sow->lambda);
	else if(H_MINUS(right[j],left[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_1]*sow->lambda);
	else
	  ZTIMES(x,&ztemp,couplings[LAMBDA_4]*sow->lambda);
      }
      if(H_EQUAL(right[j],left[j]) && H_MINUS(left[j+1],left[j+2])){
	sow=&fun13[r[j]][l[j+2]][l[j+1]];
	ZTIMES(x,&ztemp,couplings[GGN]*sow->ggn);
	if(H_EQUAL(right[j],left[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_2]*sow->lambda);
	else if(H_MINUS(right[j],left[j+1]))
	  ZTIMES(x,&ztemp,couplings[LAMBDA_1]*sow->lambda);
	else
	  ZTIMES(x,&ztemp,couplings[LAMBDA_4]*sow->lambda);
      }
      if(H_EQUAL(right[j],left[j+1]) && H_MINUS(left[j],left[j+2]) && 
	 H_ORTHO(left[j],right[j])){
	sow=&fun13[r[j]][l[j]][l[j+1]];
	ZTIMES(x,&ztemp,-couplings[BETA]*sow->lambda);
      } 
    } 
  }
  return;
}


/*  This is the interaction for the 3+1 dimensional case, including 
    a phase factor. makephaselist makes a list of angles beforehand.
    makephaselist(0,kt) should give results identical with the function
    interact_glue(...) above. The list is long enough to accomodate
    winding number up to NPMAX and 2*K upto KKMAX.

    interact_gluec(...) is NOT valid for nonzero winding number.
*/

void interact_glueballc(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings){
  int i,j;
  int l1,l2,l3,r1,r2,r3;
  partype ka,kb;
  complex_element phase,t[HINDEX];
  fun22_forms *row;
  fun13_forms *sow;
  int hl1[HINDEX],hl2[HINDEX],hl3[HINDEX],hr1[HINDEX],hr2[HINDEX],hr3[HINDEX];
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
  int il,ir;
  memcpy(leftleft,left,npl*sizeof(partype));
  memcpy(leftleft+npl,left,npl*sizeof(partype));
  memcpy(rightright,right,npr*sizeof(partype));
  memcpy(rightright+npr,right,npr*sizeof(partype));
  x->re=0.0; x->im=0.0;
  if(npl==npr && npl >1){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
	  l1=gluedecode(leftleft[il],hl1); l2=gluedecode(leftleft[il+1],hl2);
	  r1=gluedecode(rightright[ir],hr1); r2=gluedecode(rightright[ir+1],hr2);
	  for(i=0; i<HINDEX; i++){
	    j = -(2*l1+1)*hl1[i]-(2*l2+1)*(hl2[i]+2*hl1[i])
	      +(2*r1+1)*hr1[i]+(2*r2+1)*(hr2[i]+2*hr1[i]);
	    j/=2;
	    t[i].re=phaselist[abs(j)][i].re;
	    t[i].im=phaselist[abs(j)][i].im;
	    if(j<0)t[i].im *= -1.0;
	  }
	  phase.re=t[0].re*t[1].re-t[0].im*t[1].im;
	  phase.im=t[0].re*t[1].im+t[0].im*t[1].re;
	  sort_k(&ka,&kb,l1,l2,r1,r2);
	  row=&fun22[l1+l2][ka][kb];
 /* pinching terms are equal to the half and 1 the lambda_1 and lambda_4 
    interactions, respectively, for npl=npr=2*/
          if(EQUAL(hl1,hr1)){
	    ZTIMES(x,&phase,couplings[MM]*row->mass+couplings[GGN]*row->ggn);
	    if(EQUAL(hl1,hl2))
	      ZTIMES(x,&phase,couplings[LAMBDA_2]*row->lambda);
	    else if(MINUS(hl1,hl2)) {
	      ZTIMES(x,&phase,couplings[LAMBDA_1]*row->lambda);
	      if(npl==2)ZTIMES(x,&phase,0.5*couplings[LAMBDA_3]*row->lambda);
	    } else { 
	      ZTIMES(x,&phase,couplings[LAMBDA_4]*row->lambda);
	      if(npl==2 && MINUS(hl1,hl2))
		ZTIMES(x,&phase,couplings[LAMBDA_5]*row->lambda);
	    }
	  }
          if(MINUS(hl1,hl2)){
            if(npl!=2)
	      ZTIMES(x,&phase,-couplings[GGN]*row->annihilate);
	    if(EQUAL(hl1,hr1)){
	      ZTIMES(x,&phase,couplings[LAMBDA_1]*row->lambda);
	      if(npl==2)ZTIMES(x,&phase,0.5*couplings[LAMBDA_3]*row->lambda);
	    }
	    else if(MINUS(hl1,hr1))
	      ZTIMES(x,&phase,couplings[LAMBDA_2]*row->lambda);
	    else {
	      ZTIMES(x,&phase,couplings[LAMBDA_4]*row->lambda);
	       if(npl==2 && MINUS(hl1,hl2))
		 ZTIMES(x,&phase,couplings[LAMBDA_5]*row->lambda);
	    } 
	  }
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hl2))
	      ZTIMES(x,&phase,-couplings[BETA]*row->lambda);
	}
      }
   } else if(npl+2==npr){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
	  l1=gluedecode(leftleft[il],hl1); r1=gluedecode(rightright[ir],hr1); 
	  r2=gluedecode(rightright[ir+1],hr2); r3=gluedecode(rightright[ir+2],hr3);
	  for(i=0; i<HINDEX; i++){
	    j = -(2*l1+1)*hl1[i]
	       +(2*r1+1)*hr1[i]+(2*r2+1)*(hr2[i]+2*hr1[i])
	       +(2*r3+1)*(hr3[i]+2*hr2[i]+2*hr1[i]);
	    j/=2;
	    t[i].re=phaselist[abs(j)][i].re;
	    t[i].im=phaselist[abs(j)][i].im;
	    if(j<0)t[i].im *= -1.0;
	  }
	  phase.re=t[0].re*t[1].re-t[0].im*t[1].im;
	  phase.im=t[0].re*t[1].im+t[0].im*t[1].re;
          if(EQUAL(hl1,hr3)){
	    sow=&fun13[l1][r1][r2];	      
	    ZTIMES(x,&phase,couplings[GGN]*sow->ggn);
            if(EQUAL(hl1,hr2))
	      ZTIMES(x,&phase,couplings[LAMBDA_2]*sow->lambda);
	    else if(MINUS(hl1,hr2))
	      ZTIMES(x,&phase,couplings[LAMBDA_1]*sow->lambda);
	    else
	      ZTIMES(x,&phase,couplings[LAMBDA_4]*sow->lambda);
	  }
          if(EQUAL(hl1,hr1)){
	    sow=&fun13[l1][r3][r2];
	    ZTIMES(x,&phase,couplings[GGN]*sow->ggn);
            if(EQUAL(hl1,hr2))
	      ZTIMES(x,&phase,couplings[LAMBDA_2]*sow->lambda);
	    else if(MINUS(hl1,hr2))
	      ZTIMES(x,&phase,couplings[LAMBDA_1]*sow->lambda);
	    else
	      ZTIMES(x,&phase,couplings[LAMBDA_4]*sow->lambda);
	  }
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hr1)){
	    sow=&fun13[l1][r1][r2];
	    ZTIMES(x,&phase,-couplings[BETA]*sow->lambda);
	  }
        }
      } 
  } else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
	  l1=gluedecode(leftleft[il],hl1); l2=gluedecode(leftleft[il+1],hl2);
	  l3=gluedecode(leftleft[il+2],hl3);r1=gluedecode(rightright[ir],hr1);
	  for(i=0; i<HINDEX; i++){
	    j = -(2*l1+1)*hl1[i]-(2*l2+1)*(hl2[i]+2*hl1[i])
	        -(2*l3+1)*(hl3[i]+2*hl2[i]+2*hl1[i])
	        +(2*r1+1)*hr1[i];
	    j/=2;
	    t[i].re=phaselist[abs(j)][i].re;
	    t[i].im=phaselist[abs(j)][i].im;
	    if(j < 0) t[i].im *= -1.0;
	  }
	  phase.re=t[0].re*t[1].re-t[0].im*t[1].im;
	  phase.im=t[0].re*t[1].im+t[0].im*t[1].re;
          if(EQUAL(hr1,hl3)){
	    sow=&fun13[r1][l1][l2];
	    ZTIMES(x,&phase,couplings[GGN]*sow->ggn);
            if(EQUAL(hr1,hl2))
	      ZTIMES(x,&phase,couplings[LAMBDA_2]*sow->lambda);
	    else if(MINUS(hr1,hl2))
	      ZTIMES(x,&phase,couplings[LAMBDA_1]*sow->lambda);
	    else
	      ZTIMES(x,&phase,couplings[LAMBDA_4]*sow->lambda);
	  }
          if(EQUAL(hr1,hl1)){
	    sow=&fun13[r1][l3][l2];
	    ZTIMES(x,&phase,couplings[GGN]*sow->ggn);
            if(EQUAL(hr1,hl2))
	      ZTIMES(x,&phase,couplings[LAMBDA_2]*sow->lambda);
	    else if(MINUS(hr1,hl2))
	      ZTIMES(x,&phase,couplings[LAMBDA_1]*sow->lambda);
	    else
	      ZTIMES(x,&phase,couplings[LAMBDA_4]*sow->lambda);
	  }
	  if(EQUAL(hr1,hl2)&&ORTHO(hl1,hr1)){
	    sow=&fun13[r1][l1][l2];
	    ZTIMES(x,&phase,-couplings[BETA]*sow->lambda);
	  }
	}
      }
  }
/* here we take the case that the matrix element is the mass */
  else if(npl==1 && npr==1){
    x[0].re=1.0/((element) gluedecode(left[0],hl1)+0.5);
  }
  return;
}

/* 
    interact_glueballcu(x,left,right,npl,npr)

    This is a version for testing purposes which uses a
    center of mass with equal weighting for all particles.
    It's main purpose is to Check the other code.
    The angle associated with the transverse momentum
    and is obtained from the global variable "phasangleu".
    This routine can be used to calculate for nonzero winding number.
*/

void interact_glueballcu(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings){
  partype ka,kb;
  fun22_forms *row;
  fun13_forms *sow;
  int i,j,il,ir;
  int l[NPMAX+2],r[NPMAX+2],hl[NPMAX+2][HINDEX],hr[NPMAX+2][HINDEX];
  complex_element phase;
  element tl[NPMAX][HINDEX],tr[NPMAX][HINDEX],center,t; 
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
  memcpy(leftleft,left,npl*sizeof(partype));
  memcpy(leftleft+npl,left,npl*sizeof(partype));
  memcpy(rightright,right,npr*sizeof(partype));
  memcpy(rightright+npr,right,npr*sizeof(partype));
  /*  The following section calculates the distance of each 
      particle from the Center of mass.  */
  for(i=0; i<npl; i++)l[i]=gluedecode(leftleft[i],hl[i]);
  for(i=0; i<npr; i++)r[i]=gluedecode(rightright[i],hr[i]);
  for(j=0; j<HINDEX; j++){
    tl[0][j]=0.0; center=0.0;
    for(i=1; i<npl; i++)
      center += tl[i][j] = tl[i-1][j]+0.5*(element)(hl[i][j]+hl[i-1][j]);
    center *= 1.0/(element) npl;
    for(i=0; i<npl; i++)tl[i][j] -= center;
    tr[0][j]=0.0; center=0.0;  
    for(i=1; i<npr; i++)
      center+=tr[i][j]=tr[i-1][j]+0.5*(element)(hr[i][j]+hr[i-1][j]);
    center *= 1.0/(element) npr;
    for(i=0; i<npr; i++)tr[i][j] -= center;
  }
  for(i=0; i<2; i++){
    l[npl+i]=l[i]; r[npr+i]=r[i];
    for(j=0; j<HINDEX; j++){
      hl[npl+i][j]=hl[i][j]; hr[npr+i][j]=hr[i][j];
    }
  }
  /*    now go to the interactions  */
  for(i=0; i<NPARAMS; i++){x[i].re=0.0; x[i].im=0.0;}
  if(npl==npr && npl >1){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
	  for(i=0, t=0.0; i <HINDEX; i++)
	    t += phaseangleu[i]*
	      (tl[il][i]-tr[ir][i]-0.5*(element) (hl[il][i]-hr[ir][i]));
	  phase.re=cos(t); phase.im=sin(t);
	  sort_k(&ka,&kb,l[il],l[il+1],r[ir],r[ir+1]);
	  row=&fun22[l[il]+l[il+1]][ka][kb];
          if(EQUAL(hl[il],hr[ir])){
	    x[0].re+=phase.re*row->mass; x[1].re+=phase.re*row->ggn;
	    x[0].im+=phase.im*row->mass; x[1].im+=phase.im*row->ggn;
            if(EQUAL(hl[il],hl[il+1]))
	      {x[4].re+=phase.re*row->lambda; x[4].im+=phase.im*row->lambda;}
	    else if(MINUS(hl[il],hl[il+1]))
	      {x[3].re+=phase.re*row->lambda; x[3].im+=phase.im*row->lambda;}
	    else
	      {x[6].re+=phase.re*row->lambda; x[6].im+=phase.im*row->lambda;}
	  } 
          if(MINUS(hl[il],hl[il+1])){
            if(npl!=2){x[1].re-=phase.re*row->annihilate; x[1].im-=phase.im*row->annihilate;}
            if(EQUAL(hl[il],hr[ir]))
	      {x[3].re+=phase.re*row->lambda; x[3].im+=phase.im*row->lambda;}
	    else if(MINUS(hl[il],hr[ir]))
	      {x[4].re+=phase.re*row->lambda; x[4].im+=phase.im*row->lambda;}
	    else
	      {x[6].re+=phase.re*row->lambda; x[6].im+=phase.im*row->lambda;}
          } 
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hl[il+1]))
            {x[2].re-=phase.re*row->lambda; x[2].im-=phase.im*row->lambda;}
        } 
      } 
 /* pinching terms are equal to the half and 1 the lambda_1 and lambda_4 
    interactions, respectively, for npl=npr=2*/
    if(npl==2){
      x[5].re=0.5*x[3].re; x[5].im=0.5*x[3].im;
      if(MINUS(hl[0],hl[1])){
	x[7].re=x[6].re; x[7].im=x[6].im;
      }
    }
  } else if(npl+2==npr){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
	  for(i=0, t=0.0; i <HINDEX; i++)
	    t += phaseangleu[i]*
	      (tl[il][i]-tr[ir][i]-0.5*(element) (hl[il][i]-hr[ir][i]));
	  phase.re=cos(t); phase.im=sin(t);
          if(EQUAL(hl[il],hr[ir+2])){
	    sow=&fun13[l[il]][r[ir]][r[ir+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}	      
	  }
          if(EQUAL(hl[il],hr[ir])){
	    sow=&fun13[l[il]][r[ir+2]][r[ir+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}	      
	  }
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hr[ir])){
	    sow=&fun13[l[il]][r[ir]][r[ir+1]];
            {x[2].re-=phase.re*sow->lambda; x[2].im-=phase.im*sow->lambda;} 
	  }
        }
      } 
  } else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
	  for(i=0, t=0.0; i <HINDEX; i++)
	    t += phaseangleu[i]*
	      (tl[il][i]-tr[ir][i]-0.5*(element) (hl[il][i]-hr[ir][i]));
	  phase.re=cos(t); phase.im=sin(t);
          if(EQUAL(hr[ir],hl[il+2])){
	    sow=&fun13[r[ir]][l[il]][l[il+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}
	  }
          if(EQUAL(hr[ir],hl[il])){
	    sow=&fun13[r[ir]][l[il+2]][l[il+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}
	  }
	  if(EQUAL(hr[ir],hl[il+1])&&ORTHO(hl[il],hr[ir])){
	    sow=&fun13[r[ir]][l[il]][
					l[il+1]];
	    {x[2].re-=phase.re*sow->lambda; x[2].im-=phase.im*sow->lambda;}
	  }
        }
      } 
  }
/* here we take the case that the matrix element is the mass */
  else if(npl==1 && npr==1){
    x[0].re=1.0/((element) l[0]+0.5);
  }
#if UDEBUG
  if(npr==npl+2)printf("      Result: %f + %f i\n",x[3].re,x[3].im);
#endif
  return;
}

/* 
    interact_glueballcv(x,left,right,npl,npr)

    This is a version which uses a
    phase convention dependent on the first particle
    in the state and is faster than interactioncu.
    This routine can be used to calculate for nonzero winding number.
*/

void interact_glueballcv(complex_element *x, partype *left, int npl, 
		    partype *right, int npr, element *couplings){
  int i,j,il,ir;
  int l[NPMAX+2],r[NPMAX+2],hl[NPMAX+2][HINDEX],hr[NPMAX+2][HINDEX];
  complex_element phase,t[HINDEX];
  partype ka,kb;
  fun22_forms *row;
  fun13_forms *sow;
  int tl[NPMAX][HINDEX],tr[NPMAX][HINDEX];
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
  memcpy(leftleft,left,npl*sizeof(partype));
  memcpy(leftleft+npl,left,npl*sizeof(partype));
  memcpy(rightright,right,npr*sizeof(partype));
  memcpy(rightright+npr,right,npr*sizeof(partype));
  /*  The following section calculates the distance of each 
      particle from the left "edge". We set tl[0] and tr[0]
      based on the fact that we only use differences 
      tl[i]-tr[j] to calculate phases.                        */
  for(i=0; i<npl; i++)l[i]=gluedecode(leftleft[i],hl[i]);
  for(i=0; i<npr; i++)r[i]=gluedecode(rightright[i],hr[i]);
  for(j=0; j<HINDEX; j++){
    tl[0][j]=0;
    for(i=1; i<npl; i++)tl[i][j]=tl[i-1][j]+hl[i][j]+hl[i-1][j];
    tr[0][j]=hr[0][j]-hl[0][j];  
    for(i=1; i<npr; i++)tr[i][j]=tr[i-1][j]+hr[i][j]+hr[i-1][j];
  }
  for(i=0; i<2; i++){
    l[npl+i]=l[i]; r[npr+i]=r[i];
    for(j=0; j<HINDEX; j++){
      hl[npl+i][j]=hl[i][j]; hr[npr+i][j]=hr[i][j];
    }
  }
  /*    now go to the interactions  */
  for(i=0; i<NPARAMS; i++){x[i].re=0.0; x[i].im=0.0;}
  if(npl==npr && npl >1){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
	  for(i=0; i<HINDEX; i++){
	    j=tl[il][i]-tr[ir][i]-hl[il][i]+hr[ir][i];
	    t[i].re=phaselistv[abs(j)][i].re;
	    t[i].im=phaselistv[abs(j)][i].im;
	    if(j<0)t[i].im *= -1;
	  }
	  phase.re=t[0].re*t[1].re-t[0].im*t[1].im;
	  phase.im=t[0].re*t[1].im+t[0].im*t[1].re;	  
	  sort_k(&ka,&kb,l[il],l[il+1],r[ir],r[ir+1]);
	  row=&fun22[l[il]+l[il+1]][ka][kb];
          if(EQUAL(hl[il],hr[ir])){
	    x[0].re+=phase.re*row->mass; x[1].re+=phase.re*row->ggn;
	    x[0].im+=phase.im*row->mass; x[1].im+=phase.im*row->ggn;
            if(EQUAL(hl[il],hl[il+1]))
	      {x[4].re+=phase.re*row->lambda; x[4].im+=phase.im*row->lambda;}
	    else if(MINUS(hl[il],hl[il+1]))
	      {x[3].re+=phase.re*row->lambda; x[3].im+=phase.im*row->lambda;}
	    else
	      {x[6].re+=phase.re*row->lambda; x[6].im+=phase.im*row->lambda;}
	  } 
          if(MINUS(hl[il],hl[il+1])){
            if(npl!=2){x[1].re-=phase.re*row->annihilate; x[1].im-=phase.im*row->annihilate;}
            if(EQUAL(hl[il],hr[ir]))
	      {x[3].re+=phase.re*row->lambda; x[3].im+=phase.im*row->lambda;}
	    else if(MINUS(hl[il],hr[ir]))
	      {x[4].re+=phase.re*row->lambda; x[4].im+=phase.im*row->lambda;}
	    else
	      {x[6].re+=phase.re*row->lambda; x[6].im+=phase.im*row->lambda;}
          } 
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hl[il+1]))
            {x[2].re-=phase.re*row->lambda; x[2].im-=phase.im*row->lambda;}
        } 
      } 
 /* pinching terms are equal to the half and 1 the lambda_1 and lambda_4 
    interactions, respectively, for npl=npr=2*/
    if(npl==2){
      x[5].re=0.5*x[3].re; x[5].im=0.5*x[3].im;
      if(MINUS(hl[0],hl[1])){
	x[7].re=x[6].re; x[7].im=x[6].im;
      }
    }
  } else if(npl+2==npr){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
	  for(i=0; i<HINDEX; i++){
	    j=tl[il][i]-tr[ir][i]-hl[il][i]+hr[ir][i];
	    t[i].re=phaselistv[abs(j)][i].re;
	    t[i].im=phaselistv[abs(j)][i].im;
	    if(j<0)t[i].im *= -1;
	  }
	  phase.re=t[0].re*t[1].re-t[0].im*t[1].im;
	  phase.im=t[0].re*t[1].im+t[0].im*t[1].re;	  
          if(EQUAL(hl[il],hr[ir+2])){
	    sow=&fun13[l[il]][r[ir]][r[ir+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}	      
	  }
          if(EQUAL(hl[il],hr[ir])){
	    sow=&fun13[l[il]][r[ir+2]][r[ir+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}	      
	  }
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hr[ir])){
	    sow=&fun13[l[il]][r[ir]][
					r[ir+1]];
            {x[2].re-=phase.re*sow->lambda; x[2].im-=phase.im*sow->lambda;} 
          }
        }
      } 
  } else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
	  for(i=0; i<HINDEX; i++){
	    j=tl[il][i]-tr[ir][i]-hl[il][i]+hr[ir][i];
	    t[i].re=phaselistv[abs(j)][i].re;
	    t[i].im=phaselistv[abs(j)][i].im;
	    if(j<0)t[i].im *= -1;
	  }
	  phase.re=t[0].re*t[1].re-t[0].im*t[1].im;
	  phase.im=t[0].re*t[1].im+t[0].im*t[1].re;	  
          if(EQUAL(hr[ir],hl[il+2])){
	    sow=&fun13[r[ir]][l[il]][l[il+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}
	  }
          if(EQUAL(hr[ir],hl[il])){
	    sow=&fun13[r[ir]][l[il+2]][l[il+1]];
            x[1].re+=phase.re*sow->ggn; x[1].im+=phase.im*sow->ggn;
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*sow->lambda; x[4].im+=phase.im*sow->lambda;}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*sow->lambda; x[3].im+=phase.im*sow->lambda;}
	    else
	      {x[6].re+=phase.re*sow->lambda; x[6].im+=phase.im*sow->lambda;}
	  }
	  if(EQUAL(hr[ir],hl[il+1])&&ORTHO(hl[il],hr[ir])){
	    sow=&fun13[r[ir]][l[il]][l[il+1]];
	    {x[2].re-=phase.re*sow->lambda; x[2].im-=phase.im*sow->lambda;}
	  }
        }
      } 
  }
/* here we take the case that the matrix element is the mass */
  else if(npl==1 && npr==1){
    x[0].re=1.0/((element) l[0]+0.5);
  }
  return;
}


/*  
    For glueballs in 2+1 dimensions, this orders the 0^++* and 2^++ 
    states according to particle number content.  This is important 
    for correct extrapolation in 
    p-truncation and K.  Here we use the fact that the first 
    (K+1)%2 (multi=1) and K (multi=0) basis states are two particles,
    for o=1.  It does not shuffle the eigenvectors.  

    In 3+1 dimensions, it does nothing.

    BvdS, January 2002: This routine is currently *not* working.
*/

#if HINDEX==1
#define NOP 4 /* number of states to measure particle content */
#define MIN2 0.5 /*minimum fraction of two particle content expected */
                 /*before any action is taken                        */
#endif

void shuffle(full_basis *s, element *vec,
	     element *val, int nval){
#if HINDEX==1  /* only do something in 2+1 dimensions */
 int inc,top;
  element temp,p[NOP];
  size_t i,j,n;
  if(iswinding(s->ht))return;
  if((s->multi==0&&s2->multi!=0)||(s->multi==0&&s2->multi!=0)){
     fprintf(stderr,__FILE__":  only 1 multiplet=0 in shuffle, exiting\n");
    return;
  }
  if(s->multi==s2->multi&&s->o==s2->o)
    n=s->length;
  else
    n=s->length+s2->length;
  if(s->multi==0){
    inc=1; top=s->kt/2;
  } else {
    inc=0; top=(s->kt+2)/4;
  }

  if(((s->o==1&&s->multi==1)||(s2->o==1&&s2->multi==1)||
      ((s->o==1||s2->o==1)&&s->multi==0))&&nval>2+inc){
    for(j=0; j<NOP && j<nval; j++)
      for(i=0, p[j]=0; i<top; i++)
        p[j]+=vec[j*n+i]*vec[j*n+i];
    
#if 1    /* Actually do shift */
    if(p[1+inc]>p[2+inc] && p[1+inc]>MIN2){
      temp=val[1+inc];
      val[1+inc]=val[2+inc];
      val[2+inc]=temp;
    }
#endif
#if 0 /* print out warning! */
    if(p[1+inc]<MIN2 && p[2+inc]<MIN2){
      fprintf(stderr,__FILE__":  J^++ has wrong 2-particle content in shuffle,"
              " K=%i/2 p<=%i\n",s->kt,s->npmax);
      for(i=0; i<NOP && i<nval; i++)
        fprintf(stderr,"     %e %e\n",val[i],p[i]);
    }
#endif
  }
#endif
  return;
}

/* string containing couplings */

char *coupling_string(int t){
  int i;
#if HINDEX==1
  static char string[NPARAMS][20]={"mu^2","g^2 N","lambda_1","lambda_2",
				   "lambda_3","tau_1","tau_2"};
  int which[NPARAMS]={MM,GGN,LAMBDA_1,LAMBDA_2,LAMBDA_3,TAU_1,TAU_2};
#elif HINDEX==2
  static char string[NPARAMS][20]={"mu^2","g^2 N","beta","lambda_1","lambda_2",
				   "lambda_3","lambda_4","lambda_5",
				   "tau_1","tau_2","kappa_S","kappa_A","mu_F",
				   "mu_0^2","mu_1^2","mu_2^2"};
  int which[NPARAMS]={MM,GGN,BETA,LAMBDA_1,LAMBDA_2,LAMBDA_3,LAMBDA_4,
		      LAMBDA_5,TAU_1,TAU_2,KAPPA_S,KAPPA_A,
		      MU_F,MMF_0,MMF_1,MMF_2};
#endif
  for(i=0; i<NPARAMS; i++)
    if(which[i]==t)return string[i];
  return NULL;
}