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

/*  In this file are routines to imput improved matrix elements
    and calculate matrix elements of the interaction in 2+1 dimensions.
    In 2+2 dimensions, then the additional file interaction2.c is needed.
    Matrix elements for the longitudinal string tension using an open 
    string are calculated in the file interaction3.c
*/


/*
    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[KKMAX];
static element phaseangleu;
static complex_e phaselistv[NPMAX];

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

#endif




/*  multi is a basis of irreducable representations of the
    lattice symmetries */


/*  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 == 1 && HSHIFT>0

/* compare spectators using memcmp nad no modulus operation:
   somewhat faster on local machines but should be much faster
   on vector machines... */


/*  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;
  int l1,l2,l3,r1,r2,r3;
  int hl1,hl2,hl3,hr1,hr2,hr3;
  complex_e phase;
  element *row;
  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].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);
	  if(hl1 != hl2 && hl1 != hr1){       
	    phase.re=phaselist[l1+l2+1].re;
	    phase.im=phaselist[l1+l2+1].im;
	    if(hl1 > 0) phase.im *= -1.0;
	  } else if(hl1==hl2){
	    phase.re=phaselist[abs(l1-r1)].re;
	    phase.im=phaselist[abs(l1-r1)].im;
	    if(l1<r1) phase.im *= -1.0;
	    if(hl1<0) phase.im *= -1.0;
	  } else {
	    phase.re=1.0; phase.im=0.0;
	  }
	  row=twoelements[l1+l2]+four*((l1+l2+1)*l1+r1);
          if(hl1==hr1){	    
	    x[0].re+=phase.re*row[0]; x[0].im+=phase.im*row[0];
	    x[1].re+=phase.re*row[1]; x[1].im+=phase.im*row[1];
	    if(hl1==hl2)
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else
	      {x[2].re+=phase.re*row[3]; x[2].im+=phase.im*row[3];}
	  } 
          if(hl1 != hl2){
            if(npl!=2)
	      {x[1].re-=phase.re*row[2]; x[1].im-=phase.im*row[2];}
            if(hl1==hr1)
	      {x[2].re+=phase.re*row[3]; x[2].im+=phase.im*row[3];}
	    else
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	  } 
        } 
      } 
    /* pinching term is equal to half the lambda_1 interaction for npl=npr=2*/
      if(npl==2){x[4].re=0.5*x[2].re; x[4].im=0.5*x[2].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);
	  if(hl1 == hr2){
	    if(hl1 == hr1){
	      phase.re=phaselist[r2+r3+1].re;
	      phase.im=phaselist[r2+r3+1].im;
	      if(hl1 < 0) phase.im *= -1.0;/* < */
	    } else {
	      phase.re=phaselist[r1+r2+1].re;
	      phase.im=phaselist[r1+r2+1].im;	      
	      if(hl1 > 0) phase.im *= -1.0;/* > */
	    }
	  } else {
	    phase.re=1.0;
	    phase.im=0.0;
	  } 
          if(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(hl1==hr2)
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[2].re+=phase.re*row[1]; x[2].im+=phase.im*row[1];}
	  }
          if(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(hl1==hr2)
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {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);
	  if(hr1 == hl2){
	    if(hr1 == hl1){
	      phase.re=phaselist[l2+l3+1].re;
	      phase.im=phaselist[l2+l3+1].im;
	      if(hr1>0) phase.im *= -1.0; /* > */
	    } else {
	      phase.re=phaselist[l1+l2+1].re;
	      phase.im=phaselist[l1+l2+1].im;	      
	      if(hr1 < 0) phase.im *= -1.0; /* < */
	    }
	  } else {
	    phase.re=1.0; phase.im=0.0;
	  } 
	  if(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(hr1==hl2)
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[2].re+=phase.re*row[1]; x[2].im+=phase.im*row[1];}
	  } 
          if(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(hr1==hl2)
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {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(leftleft[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.
*/

#define UDEBUG 1==0

void interactioncu(complex_e *x, 
		 partype *left, partype *right, int npl, int npr){
  int i,il,ir;
  int l[NPMAX+2],r[NPMAX+2],hl[NPMAX+2],hr[NPMAX+2];
  complex_e phase;
  element tl[NPMAX],tr[NPMAX],center,*row;
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
#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));
  /*  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);
  tl[0]=0.0;
  for(i=1, center=0; i<npl; i++)
    center+=tl[i]=tl[i-1]+((element)(hl[i]+hl[i-1]))/2.0;
  center *= 1.0/(element) npl;
  for(i=0; i<npl; i++)tl[i]-=center;
  for(i=0; i<npr; i++)r[i]=hitol(rightright[i],hr+i);
  tr[0]=0.0;  
  for(i=1, center=0; i<npr; i++)
    center+=tr[i]=tr[i-1]+((element)(hr[i]+hr[i-1]))/2.0;
  center *= 1.0/(element) npr;
  for(i=0; i<npr; i++)tr[i]-=center;
  for(i=0; i<2; i++){
    l[npl+i]=l[i];  hl[npl+i]=hl[i];
    r[npr+i]=r[i];  hr[npr+i]=hr[i];
  }
  /*    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
	  center=tl[il]-tr[ir]-0.5*(element) (hl[il]-hr[ir]);
	  phase.re=cos(phaseangleu*center);
	  phase.im=sin(phaseangleu*center);
	  row=twoelements[l[il]+l[il+1]]+four*((l[il]+l[il+1]+1)*l[il]+r[ir]);
          if(hl[il]==hr[ir]){	    
	    x[0].re+=phase.re*row[0]; x[0].im+=phase.im*row[0];
	    x[1].re+=phase.re*row[1]; x[1].im+=phase.im*row[1];
	    if(hl[il]==hl[il+1])
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else
	      {x[2].re+=phase.re*row[3]; x[2].im+=phase.im*row[3];}
	  } 
          if(hl[il] != hl[il+1]){
            if(npl!=2)
	      {x[1].re-=phase.re*row[2]; x[1].im-=phase.im*row[2];}
            if(hl[il]==hr[ir])
	      {x[2].re+=phase.re*row[3]; x[2].im+=phase.im*row[3];}
	    else
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	  } 
        } 
      }
    /* pinching term is equal to half the lambda_1 interaction for npl=npr=2*/
    if(npl==2){x[4].re=0.5*x[2].re; x[4].im=0.5*x[2].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
	  center=tl[il]-tr[ir]-0.5*(element) (hl[il]-hr[ir]);
	  phase.re=cos(phaseangleu*center);
	  phase.im=sin(phaseangleu*center);
          if(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(hl[il]==hr[ir+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[2].re+=phase.re*row[1]; x[2].im+=phase.im*row[1];}
	  }
          if(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(hl[il]==hr[ir+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {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
	  center=tl[il]-tr[ir]-0.5*(element) (hl[il]-hr[ir]);
	  phase.re=cos(phaseangleu*center);
	  phase.im=sin(phaseangleu*center);
	  if(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(hr[ir]==hl[il+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[2].re+=phase.re*row[1]; x[2].im+=phase.im*row[1];}
	  } 
          if(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(hr[ir]==hl[il+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {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 it 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,il,ir;
  int l[NPMAX+2],r[NPMAX+2],hl[NPMAX+2],hr[NPMAX+2];
  complex_e phase;
  element *row;
  int tl[NPMAX],tr[NPMAX];
  partype leftleft[NPMAX*2],rightright[NPMAX*2];
#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));
  /*  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);
  tl[0]=0;
  for(i=1; i<npl; i++)tl[i]=tl[i-1]+(hl[i]+hl[i-1])/2;
  for(i=0; i<npr; i++)r[i]=hitol(rightright[i],hr+i);
  tr[0]=(hr[0]-hl[0])/2;  
  for(i=1; i<npr; i++)tr[i]=tr[i-1]+(hr[i]+hr[i-1])/2;
  for(i=0; i<2; i++){
    l[npl+i]=l[i];  hl[npl+i]=hl[i];
    r[npr+i]=r[i];  hr[npr+i]=hr[i];
  }
  /*    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
	  i=tl[il]-tr[ir]-(hl[il]-hr[ir])/2;
	  phase.re=phaselistv[abs(i)].re;
	  phase.im=phaselistv[abs(i)].im;
	  if(i<0)phase.im *= -1;
	  row=twoelements[l[il]+l[il+1]]+four*((l[il]+l[il+1]+1)*l[il]+r[ir]);
          if(hl[il]==hr[ir]){	    
	    x[0].re+=phase.re*row[0]; x[0].im+=phase.im*row[0];
	    x[1].re+=phase.re*row[1]; x[1].im+=phase.im*row[1];
	    if(hl[il]==hl[il+1])
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	    else
	      {x[2].re+=phase.re*row[3]; x[2].im+=phase.im*row[3];}
	  } 
          if(hl[il] != hl[il+1]){
            if(npl!=2)
	      {x[1].re-=phase.re*row[2]; x[1].im-=phase.im*row[2];}
            if(hl[il]==hr[ir])
	      {x[2].re+=phase.re*row[3]; x[2].im+=phase.im*row[3];}
	    else
	      {x[3].re+=phase.re*row[3]; x[3].im+=phase.im*row[3];}
	  } 
        } 
      } 
    /* pinching term is equal to half the lambda_1 interaction for npl=npr=2*/
    if(npl==2){x[4].re=0.5*x[2].re; x[4].im=0.5*x[2].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
	  i=tl[il]-tr[ir]-(hl[il]-hr[ir])/2;
	  phase.re=phaselistv[abs(i)].re;
	  phase.im=phaselistv[abs(i)].im;
	  if(i<0)phase.im *= -1;
          if(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(hl[il]==hr[ir+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[2].re+=phase.re*row[1]; x[2].im+=phase.im*row[1];}
	  }
          if(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(hl[il]==hr[ir+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {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
	  i=tl[il]-tr[ir]-(hl[il]-hr[ir])/2;
	  phase.re=phaselistv[abs(i)].re;
	  phase.im=phaselistv[abs(i)].im;
	  if(i<0)phase.im *= -1;
	  if(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(hr[ir]==hl[il+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {x[2].re+=phase.re*row[1]; x[2].im+=phase.im*row[1];}
	  } 
          if(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(hr[ir]==hl[il+1])
	      {x[3].re+=phase.re*row[1]; x[3].im+=phase.im*row[1];}
	    else
	      {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 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.  */

#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                        */
 
void shuffle(sectors *s, sectors *s2,element *vec,
	     element *val,int nval){
  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,"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];
    /*Actually do shift  */
#if 1==1
    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 1==0 /* print out warning! */
    if(p[1+inc]<MIN2 && p[2+inc]<MIN2){
      printf("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++)
	printf("     %e %e\n",val[i],p[i]);
    }
#endif
  }
  return;
}

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

void interaction(element *x, partype *left, partype *right, int npl, int npr){
  int i;
  int l1,l2,l3,r1,r2,r3,hl1,hl2,hl3,hr1,hr2,hr3;
  element *row;
  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(hl1 == hr1){
            x[0]+=row[0]; 
	    x[1]+=row[1];
            if(hl1==hl2)
	      x[3]+=row[3];
	    else
	      x[2]+=row[3];
	  } 
          if(hl1 != hl2){
            if(npl!=2)
	      x[1]-=row[2]; 
            if(hl1==hr1)
	      x[2]+=row[3]; 
	    else
	      x[3]+=row[3]; 
          } 
        }
      } 
    /* pinching term is equal to half the lambda_1 interaction for npl=npr=2*/
    if(npl==2)x[4]=0.5*x[2];
  } 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(hl1 == hr3){
	    row=threeelements[l1-1]+((((2*l1-r1+1)*r1)>>1)+r2)*two;
            x[1]+=row[0];
            if(hl1==hr2)
	      x[3]+=row[1];
	    else
	      x[2]+=row[1];
	  }
          if(hl1 == hr1){
	    row=threeelements[l1-1]+((((2*l1-r3+1)*r3)>>1)+r2)*two;
            x[1]+=row[0];
            if(hl1==hr2)
	      x[3]+=row[1];
	    else
	      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(hr1 == hl3){
	    row=threeelements[r1-1]+((((2*r1-l1+1)*l1)>>1)+l2)*two;
            x[1]+=row[0];
            if(hr1==hl2)
	      x[3]+=row[1];
	    else
	      x[2]+=row[1];
	  }
          if(hr1 == hl1){
	    row=threeelements[r1-1]+((((2*r1-l3+1)*l3)>>1)+l2)*two;
            x[1]+=row[0];
            if(hr1==hl2)
	      x[3]+=row[1];
	    else
	      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;
}

#endif