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


/*
  Below are subroutines for the 3+1 case only.

  multi is a basis of irreducable representations of the
  lattice symmetries 

  Define lists of phases used by various subroutines
  Here, kt = 2*K.  Half integers are defined for
  phaselist and phaselistv when HINDEX>1.
*/

#if HINDEX > 1

static complex_e phaselist[2*KKMAX][HINDEX];
static element phaseangleu[HINDEX];
static complex_e 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
}
#endif



/*  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.

*/

 
#if HINDEX == 2 && HSHIFT>0

/*  This is the interaction for the 3+1 dimensional case.
    It returns values in the vector x. */


void interaction(element *x, partype *left, partype *right, int npl, int npr){
  int i,l1,l2,l3,r1,r2,r3;
  element *row;
  int hl1[HINDEX],hl2[HINDEX],hl3[HINDEX],hr1[HINDEX],hr2[HINDEX],hr3[HINDEX];
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
  int il,ir;
#if !USEMEMCMP
  int j;
#endif
  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));
  for(i=0; i<NPARAMS; i++)x[i]=0.0;
  if(npl==npr && npl >1){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+2,rightright+ir+2,(npl-2)*sizeof(partype))){
#else
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
#endif
	  l1=hitol(leftleft[il],hl1); l2=hitol(leftleft[il+1],hl2);
	  r1=hitol(rightright[ir],hr1); r2=hitol(rightright[ir+1],hr2);
	  row=twoelements[l1+l2]+four*((l1+l2+1)*l1+r1);
          if(EQUAL(hl1,hr1)){
	    x[0]+=row[0]; x[1]+=row[1];
            if(EQUAL(hl1,hl2))
	      x[4]+=row[3];
	    else if(MINUS(hl1,hl2))
	      x[3]+=row[3];
	    else
	      x[6]+=row[3];
	  } 
          if(MINUS(hl1,hl2)){
            if(npl!=2)
	      x[1]-=row[2];
            if(EQUAL(hl1,hr1))
	      x[3]+=row[3];
	    else if(MINUS(hl1,hr1))
	      x[4]+=row[3];
	    else
	      x[6]+=row[3];
          } 
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hl2))
	    x[2]-=row[3];
        } 
      } 
 /* 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]=0.5*x[3];
      /* This is to ensure that lambda_5 does not 
         act on winding modes. */
      l1=hitol(leftleft[0],hl1); l2=hitol(leftleft[1],hl2);
      if(MINUS(hl1,hl2))
	x[7]=x[6];
    }
  } else if(npl+2==npr){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+1,rightright+ir+3,(npl-1)*sizeof(partype))){
#else
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
#endif
	  l1=hitol(leftleft[il],hl1); r1=hitol(rightright[ir],hr1); 
	  r2=hitol(rightright[ir+1],hr2); r3=hitol(rightright[ir+2],hr3);
          if(EQUAL(hl1,hr3)){
	    row=threeelements[l1-1]+((((2*l1-r1+1)*r1)>>1)+r2)*two;
	    x[1]+=row[0];
            if(EQUAL(hl1,hr2))
	      x[4]+=row[1];
	    else if(MINUS(hl1,hr2))
	      x[3]+=row[1];
	    else
	      x[6]+=row[1];
	  }
          if(EQUAL(hl1,hr1)){
	    row=threeelements[l1-1]+((((2*l1-r3+1)*r3)>>1)+r2)*two;
	    x[1]+=row[0];
            if(EQUAL(hl1,hr2))
	      x[4]+=row[1];
	    else if(MINUS(hl1,hr2))
	      x[3]+=row[1];
	    else
	      x[6]+=row[1];
	  }
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hr1)){
	    row=threeelements[l1-1]+((((2*l1-r1+1)*r1)>>1)+r2)*two;
	    x[2]-=row[1];
	  }
	} 
      } 
  }
  else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+3,rightright+ir+1,(npr-1)*sizeof(partype))){
#else
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
#endif
	  l1=hitol(leftleft[il],hl1); l2=hitol(leftleft[il+1],hl2);
	  l3=hitol(leftleft[il+2],hl3);r1=hitol(rightright[ir],hr1);
          if(EQUAL(hr1,hl3)){
	    row=threeelements[r1-1]+((((2*r1-l1+1)*l1)>>1)+l2)*two;
	    x[1]+=row[0];
            if(EQUAL(hr1,hl2))
	      x[4]+=row[1];
	    else if(MINUS(hr1,hl2))
	      x[3]+=row[1];
	    else
	      x[6]+=row[1];
	  }
          if(EQUAL(hr1,hl1)){
	    row=threeelements[r1-1]+((((2*r1-l3+1)*l3)>>1)+l2)*two;
	    x[1]+=row[0];
            if(EQUAL(hr1,hl2))
	      x[4]+=row[1];
	    else if(MINUS(hr1,hl2))
	      x[3]+=row[1];
	    else
	      x[6]+=row[1];
	  }
	  if(EQUAL(hr1,hl2)&&ORTHO(hl1,hr1)){
	    row=threeelements[r1-1]+((((2*r1-l1+1)*l1)>>1)+l2)*two;
	    x[2]-=row[1];
	  }
	} 
      } 
  }
/* here we take the case that the matrix element is the mass */
  else if(npl==1 && npr==1){
    x[0]=1.0/((element) hitol(left[0],hl1)+0.5);
  }
  return;
}

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

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

void interactionc(complex_e *x, 
		 partype *left, partype *right, int npl, int npr){
  int i,j;
  int l1,l2,l3,r1,r2,r3;
  complex_e phase,t[HINDEX];
  element *row;
  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));
  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++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+2,rightright+ir+2,(npl-2)*sizeof(partype))){
#else
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
#endif
	  l1=hitol(leftleft[il],hl1); l2=hitol(leftleft[il+1],hl2);
	  r1=hitol(rightright[ir],hr1); r2=hitol(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;
	  row=twoelements[l1+l2]+four*((l1+l2+1)*l1+r1);
          if(EQUAL(hl1,hr1)){
	    x[0].re+=phase.re*row[0]; x[1].re+=phase.re*row[1];
	    x[0].im+=phase.im*row[0]; x[1].im+=phase.im*row[1];
            if(EQUAL(hl1,hl2))
	      {x[4].re+=phase.re*row[3]; x[4].im+=phase.im*row[3];}
	    else if(MINUS(hl1,hl2))
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else
	      {x[6].re+=phase.re*row[3]; x[6].im+=phase.im*row[3];}
	  } 
          if(MINUS(hl1,hl2)){
            if(npl!=2)
	      {x[1].re-=phase.re*row[2]; x[1].im-=phase.im*row[2];}
            if(EQUAL(hl1,hr1))
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else if(MINUS(hl1,hr1))
	      {x[4].re+=phase.re*row[3]; x[4].im+=phase.im*row[3];}
	    else
	      {x[6].re+=phase.re*row[3]; x[6].im+=phase.im*row[3];}
          } 
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hl2))
            {x[2].re-=phase.re*row[3]; x[2].im-=phase.im*row[3];}
	}
      }
 /* 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;
      l1=hitol(leftleft[0],hl1); l2=hitol(leftleft[1],hl2);
      if(MINUS(hl1,hl2)){
	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++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+1,rightright+ir+3,(npl-1)*sizeof(partype))){
#else
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
#endif
	  l1=hitol(leftleft[il],hl1); r1=hitol(rightright[ir],hr1); 
	  r2=hitol(rightright[ir+1],hr2); r3=hitol(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)){
	    row=threeelements[l1-1]+((((2*l1-r1+1)*r1)>>1)+r2)*two;
	    x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hl1,hr2))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hl1,hr2))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
          if(EQUAL(hl1,hr1)){
	    row=threeelements[l1-1]+((((2*l1-r3+1)*r3)>>1)+r2)*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hl1,hr2))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hl1,hr2))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}	      
	  }
          if(EQUAL(hl1,hr2)&&ORTHO(hl1,hr1)){
	    row=threeelements[l1-1]+((((2*l1-r1+1)*r1)>>1)+r2)*two;
	    x[2].re-=phase.re*row[1]; x[2].im-=phase.im*row[1];	      
	  }
        }
      } 
  } else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+3,rightright+ir+1,(npr-1)*sizeof(partype))){
#else
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
#endif
	  l1=hitol(leftleft[il],hl1); l2=hitol(leftleft[il+1],hl2);
	  l3=hitol(leftleft[il+2],hl3);r1=hitol(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)){
	    row=threeelements[r1-1]+((((2*r1-l1+1)*l1)>>1)+l2)*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hr1,hl2))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hr1,hl2))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
          if(EQUAL(hr1,hl1)){
	    row=threeelements[r1-1]+((((2*r1-l3+1)*l3)>>1)+l2)*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hr1,hl2))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hr1,hl2))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
	  if(EQUAL(hr1,hl2)&&ORTHO(hl1,hr1)){
	    row=threeelements[r1-1]+((((2*r1-l1+1)*l1)>>1)+l2)*two;
	    {x[2].re-=phase.re*row[1]; x[2].im-=phase.im*row[1];}
	  }
	}
      }
  }
/* here we take the case that the matrix element is the mass */
  else if(npl==1 && npr==1){
    x[0].re=1.0/((element) hitol(left[0],hl1)+0.5);
  }
  return;
}

/* 
    interactioncu(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 interactioncu(complex_e *x, 
		 partype *left, partype *right, int npl, int npr){
  int i,j,il,ir;
  int l[NPMAX+2],r[NPMAX+2],hl[NPMAX+2][HINDEX],hr[NPMAX+2][HINDEX];
  complex_e phase;
  element tl[NPMAX][HINDEX],tr[NPMAX][HINDEX],center,*row,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]=hitol(leftleft[i],hl[i]);
  for(i=0; i<npr; i++)r[i]=hitol(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++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+2,rightright+ir+2,(npl-2)*sizeof(partype))){
#else
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
#endif
	  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);
	  row=twoelements[l[il]+l[il+1]]+four*((l[il]+l[il+1]+1)*l[il]+r[ir]);
          if(EQUAL(hl[il],hr[ir])){
	    x[0].re+=phase.re*row[0]; x[1].re+=phase.re*row[1];
	    x[0].im+=phase.im*row[0]; x[1].im+=phase.im*row[1];
            if(EQUAL(hl[il],hl[il+1]))
	      {x[4].re+=phase.re*row[3]; x[4].im+=phase.im*row[3];}
	    else if(MINUS(hl[il],hl[il+1]))
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else
	      {x[6].re+=phase.re*row[3]; x[6].im+=phase.im*row[3];}
	  } 
          if(MINUS(hl[il],hl[il+1])){
            if(npl!=2){x[1].re-=phase.re*row[2]; x[1].im-=phase.im*row[2];}
            if(EQUAL(hl[il],hr[ir]))
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else if(MINUS(hl[il],hr[ir]))
	      {x[4].re+=phase.re*row[3]; x[4].im+=phase.im*row[3];}
	    else
	      {x[6].re+=phase.re*row[3]; x[6].im+=phase.im*row[3];}
          } 
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hl[il+1]))
            {x[2].re-=phase.re*row[3]; x[2].im-=phase.im*row[3];}
        } 
      } 
 /* 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++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+1,rightright+ir+3,(npl-1)*sizeof(partype))){
#else
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
#endif
	  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])){
	    row=threeelements[l[il]-1]+((((2*l[il]-r[ir]+1)*r[ir])>>1)+
					r[ir+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}	      
	  }
          if(EQUAL(hl[il],hr[ir])){
	    row=threeelements[l[il]-1]+((((2*l[il]-r[ir+2]+1)*r[ir+2])>>1)+
					r[ir+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}	      
	  }
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hr[ir])){
	    row=threeelements[l[il]-1]+((((2*l[il]-r[ir]+1)*r[ir])>>1)+r[ir+1])*two;
            {x[2].re-=phase.re*row[1]; x[2].im-=phase.im*row[1];} 
	  }
        }
      } 
  } else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+3,rightright+ir+1,(npr-1)*sizeof(partype))){
#else
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
#endif
	  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])){
	    row=threeelements[r[ir]-1]+((((2*r[ir]-l[il]+1)*l[il])>>1)+
					l[il+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
          if(EQUAL(hr[ir],hl[il])){
	    row=threeelements[r[ir]-1]+((((2*r[ir]-l[il+2]+1)*l[il+2])>>1)+
					l[il+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
	  if(EQUAL(hr[ir],hl[il+1])&&ORTHO(hl[il],hr[ir])){
	    row=threeelements[r[ir]-1]+((((2*r[ir]-l[il]+1)*l[il])>>1)+
					l[il+1])*two;
	    {x[2].re-=phase.re*row[1]; x[2].im-=phase.im*row[1];}
	  }
        }
      } 
  }
/* 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;
}

/* 
    interactioncv(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 interactioncv(complex_e *x, 
		 partype *left, partype *right, int npl, int npr){
  int i,j,il,ir;
  int l[NPMAX+2],r[NPMAX+2],hl[NPMAX+2][HINDEX],hr[NPMAX+2][HINDEX];
  complex_e phase,t[HINDEX];
  element *row;
  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]=hitol(leftleft[i],hl[i]);
  for(i=0; i<npr; i++)r[i]=hitol(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++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+2,rightright+ir+2,(npl-2)*sizeof(partype))){
#else
	for(j=2; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j])break;
        if(j==npl){
#endif
	  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;	  
	  row=twoelements[l[il]+l[il+1]]+four*((l[il]+l[il+1]+1)*l[il]+r[ir]);
          if(EQUAL(hl[il],hr[ir])){
	    x[0].re+=phase.re*row[0]; x[1].re+=phase.re*row[1];
	    x[0].im+=phase.im*row[0]; x[1].im+=phase.im*row[1];
            if(EQUAL(hl[il],hl[il+1]))
	      {x[4].re+=phase.re*row[3]; x[4].im+=phase.im*row[3];}
	    else if(MINUS(hl[il],hl[il+1]))
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else
	      {x[6].re+=phase.re*row[3]; x[6].im+=phase.im*row[3];}
	  } 
          if(MINUS(hl[il],hl[il+1])){
            if(npl!=2){x[1].re-=phase.re*row[2]; x[1].im-=phase.im*row[2];}
            if(EQUAL(hl[il],hr[ir]))
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else if(MINUS(hl[il],hr[ir]))
	      {x[4].re+=phase.re*row[3]; x[4].im+=phase.im*row[3];}
	    else
	      {x[6].re+=phase.re*row[3]; x[6].im+=phase.im*row[3];}
          } 
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hl[il+1]))
            {x[2].re-=phase.re*row[3]; x[2].im-=phase.im*row[3];}
        } 
      } 
 /* 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++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+1,rightright+ir+3,(npl-1)*sizeof(partype))){
#else
	for(j=1; j<npl; j++)
	  if(leftleft[il+j] != rightright[ir+j+2])break;
        if(j==npl){
#endif
	  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])){
	    row=threeelements[l[il]-1]+((((2*l[il]-r[ir]+1)*r[ir])>>1)+
					r[ir+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}	      
	  }
          if(EQUAL(hl[il],hr[ir])){
	    row=threeelements[l[il]-1]+((((2*l[il]-r[ir+2]+1)*r[ir+2])>>1)+
					r[ir+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hl[il],hr[ir+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hl[il],hr[ir+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}	      
	  }
          if(EQUAL(hl[il],hr[ir+1])&&ORTHO(hl[il],hr[ir])){
	    row=threeelements[l[il]-1]+((((2*l[il]-r[ir]+1)*r[ir])>>1)+
					r[ir+1])*two;
            {x[2].re-=phase.re*row[1]; x[2].im-=phase.im*row[1];} 
          }
        }
      } 
  } else if(npr+2==npl){
    for(il=0; il<npl; il++)
      for(ir=0; ir<npr; ir++){
#if USEMEMCMP
        if(!memcmp(leftleft+il+3,rightright+ir+1,(npr-1)*sizeof(partype))){
#else
	for(j=1; j<npr; j++)
	  if(leftleft[il+j+2] != rightright[ir+j])break;
        if(j==npr){
#endif
	  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])){
	    row=threeelements[r[ir]-1]+((((2*r[ir]-l[il]+1)*l[il])>>1)+
					l[il+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
          if(EQUAL(hr[ir],hl[il])){
	    row=threeelements[r[ir]-1]+((((2*r[ir]-l[il+2]+1)*l[il+2])>>1)+
					l[il+1])*two;
            x[1].re+=phase.re*row[0]; x[1].im+=phase.im*row[0];
            if(EQUAL(hr[ir],hl[il+1]))
	      {x[4].re+=phase.re*row[1]; x[4].im+=phase.im*row[1];}
	    else if(MINUS(hr[ir],hl[il+1]))
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[6].re+=phase.re*row[1]; x[6].im+=phase.im*row[1];}
	  }
	  if(EQUAL(hr[ir],hl[il+1])&&ORTHO(hl[il],hr[ir])){
	    row=threeelements[r[ir]-1]+((((2*r[ir]-l[il]+1)*l[il])>>1)+
					l[il+1])*two;
	    {x[2].re-=phase.re*row[1]; x[2].im-=phase.im*row[1];}
	  }
        }
      } 
  }
/* 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;
}


/*  This currently does nothing     */

#define NOP 4 /* number of states to measure particle content */

void shuffle(sectors *s, sectors *s2,element *vec,
	     element *val,int nval){
  return;
}


#endif