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

#define PI 3.1415926535897932385

#ifndef YUKAWA_FOUR_POINT
#define YUKAWA_FOUR_POINT 0
#endif

/*    
      Matrix elements of P^- for mesons

      The quark interactions are calculated "in line"  There
      are several reasons for this:

       1. This allows quickest development based on Mma code

       2. I have not developed an "improvement" scheme, which is 
       the motivation for calculating the interactions separately.

       3. Storage of the interaction would take up a significant
       amount of memory, due to the larger values of K involved and the
       large number of functional forms of the various interactions.


To convert the Mathematica expressions, I used
emacs Query replace regexp: (\([^()]+\))\^2   with          pow(\1,2)
                            Sqrt\[\([^[]+\)\]     with      sqrt(\1)
                            \([0-9a-zA-A\])]\) \([a-zA-Z]\) \1*\2
                                  
*/


/*
  The following table was created with Mathematica using:

    <<makeham.m
    expr1=Table[Map[Table[zetaf[b[0,h1],a[0,#1],d[0,h2]],
    {h2,1/2,-1/2,-1}]&,{1,-1,2,-2}],{h1,1/2,-1/2,-1}]
    InputForm[Map[{Re[#],Im[#]}&,expr1,{3}]]

*/

const complex_element zeta[2][4][2]=
{{{{0, 0}, {0, -1}}, {{0, 0}, {0, 1}}, {{0, 0}, {-1, 0}}, {{0, 0}, {1, 0}}}, 
 {{{0, 1}, {0, 0}}, {{0, -1}, {0, 0}}, {{-1, 0}, {0, 0}}, {{1, 0}, {0, 0}}}};

#define zeta(A,B,C) zeta[(A)>>HSHIFT][(B)>>HSHIFT][(C)>>HSHIFT]


/*
  The following table was created with Mathematica using:

    <<makeham.m
    expr2=Table[Outer[Table[zetafzetaf[b[0,h1],a[0,#1],a[0,#2],d[0,h2]],
    {h2,1/2,-1/2,-1}]&,{1,-1,2,-2},{1,-1,2,-2}],{h1,1/2,-1/2,-1}] 
    InputForm[Map[{Re[#],Im[#]}&,expr2,{4}]]

*/

complex_element zetazeta[2][4][4][2]=
{{{{{1, 0}, {0, 0}}, {{-1, 0}, {0, 0}}, {{0, 1}, {0, 0}}, {{0, -1}, {0, 0}}}, 
  {{{-1, 0}, {0, 0}}, {{1, 0}, {0, 0}}, {{0, -1}, {0, 0}}, {{0, 1}, {0, 0}}}, 
  {{{0, -1}, {0, 0}}, {{0, 1}, {0, 0}}, {{1, 0}, {0, 0}}, {{-1, 0}, {0, 0}}}, 
  {{{0, 1}, {0, 0}}, {{0, -1}, {0, 0}}, {{-1, 0}, {0, 0}}, 
   {{1, 0}, {0, 0}}}}, {{{{0, 0}, {1, 0}}, {{0, 0}, {-1, 0}}, 
   {{0, 0}, {0, -1}}, {{0, 0}, {0, 1}}}, {{{0, 0}, {-1, 0}}, 
   {{0, 0}, {1, 0}}, {{0, 0}, {0, 1}}, {{0, 0}, {0, -1}}}, 
  {{{0, 0}, {0, 1}}, {{0, 0}, {0, -1}}, {{0, 0}, {1, 0}}, {{0, 0}, {-1, 0}}}, 
  {{{0, 0}, {0, -1}}, {{0, 0}, {0, 1}}, {{0, 0}, {-1, 0}}, {{0, 0}, {1, 0}}}}};

#define zetazeta(A,B,C,D) zetazeta[(A)>>HSHIFT][(B)>>HSHIFT]\
                          [(C)>>HSHIFT][D>>HSHIFT]

element elfkappa(partype k){
  element sum=0;
  int i;
  k&=KMASK;
#if 1  /* My way, method for summer 2002 */
  for(i=0; i<k; i++)sum+=1.0/QE(i);
#else /* Simon's way, from summer 2001 */
  if(k>0)
    for(i=0; i<k; i++)sum+=1.0/((element) i+1);
#endif
  return sum/(2.0*PI*QE(k));
}

element elfquarkquark(partype k1, partype k2){
  element sum=0;
  int i;
  k1&=KMASK; k2&=KMASK;
  for(i=0; i<=k1+k2; i++)
    if(i!=k1)
      sum+=1.0/pow(QE(k1)-QE(i),2);
  return sum/(4.0*PI);
}

/*  elfquarkglue[k1,k2]:  k1 is the quark {1/2,3/2,...} and k2 is the link
field {1,2,3,...} */

element elfquarkglue(partype k1, partype k2){
  element sum;
  int i;
  k1&=KMASK; k2&=KMASK;
  sum=2.0*(sqrt(1+QE(k1)/GE(k2))-1.0)/QE(k1);
  /* the upper limit depends on link zer modes */
  for(i=0; i<k1+k2; i++)
    if(i!=k1)
      sum+=(QE(k1)+2.0*GE(k2)-QE(i))/
	(sqrt(GE(k2)*(QE(k1)+GE(k2)-QE(i)))*pow(QE(k1)-QE(i),2));
  return sum/(8.0*PI);
}

/*
  find phase factor for given transition, with quark at beginning
*/

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

/* 
  This calculates matrix elements of the Hamiltonian for a meson.

  It does either zero or nonzero transverse momentum
  based on the value of current->momflag.
*/

void interact_meson(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings){
  partype left_flip[NPMAX],right_flip[NPMAX];

  x->re=0.0; x->im=0.0;

#if 1      /* for debugging only -- turn off link-link interactions */
  if(left[0]==right[0] && left[npl-1]==right[npr-1]){ /* check quarks */
    if(working.momflag)
      interact_gluec(x,left+1,npl-2,right+1,npr-2,couplings);
    else
      interact_glue(&(x->re),left+1,npl-2,right+1,npr-2,couplings);
  }
#else
  fprintf(stderr,__FILE__":  link-link interactions turned off.\n"); 
#endif 

  /*  In the following, we calculate the contibutions from the ends of
      the strings.   */
  
#if 1   /* debugging flag to turn off ends of the strings */

  /* quark-antiquark gauge interactions */
  if(npl==2 && npr==2 && H_EQUAL(left[0],right[0]) && 
     H_EQUAL(left[1],right[1])){
    if(left[0]!=right[0]){
      x->re-=couplings[GGN]/(4*PI*pow((QE(left[0])-QE(right[0])),2));
    } else {
      x->re+=couplings[GGN]*elfquarkquark(left[0],left[1]);
    }
  }

  if(working.momflag)
    interact_quarkc(x,left,npl,right,npr,couplings,npl-2);
  else
    interact_quark(x,left,npl,right,npr,couplings,npl-2);

  /* Do the anti-quark interactions by transforming the state */
  charge_conjugation(left_flip,left,TYPE_MESON,npl);
  charge_conjugation(right_flip,right,TYPE_MESON,npr);
  if(working.momflag)
    interact_quarkc(x,left_flip,npl,right_flip,npr,couplings,npl-2);
  else
    interact_quark(x,left_flip,npl,right_flip,npr,couplings,npl-2);

#endif  /* end of debugging flag to turn off end of string interactions */

#if 1 /* debug print */
  if(x->re==0.0 || fabs(x->re)<1.0e10) return;
  printf("   last x=%f%+f*I:  ",x->re,x->im);
  printstate(left,TYPE_MESON,npl);
  printf(" - ");
  printstate(right,TYPE_MESON,npr);
  printf("\n");
#endif  

  return;
}

/*
   Interactions with quarks and gluons (one end of the string)
   quark-antiquark interactions not included.

   x is not initialized

   ngl is number of gluons in left state.
*/

void interact_quark(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings, int ngl){
  int j;
  element temp;

  /* Mass terms and Yukawa self-intertias */

  if(npl==npr){
    for(j=0; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      temp=pow(couplings[MU_F],2);
      if(npl==working.nptrunc)
	temp+=couplings[MMF_0];
      if(npl==working.nptrunc-1)
	temp+=couplings[MMF_1];
      if(npl<=working.nptrunc-2)
	temp+=couplings[MMF_2];
       x->re+=temp/(2.0*QE(left[0]));
      if(npl<working.nptrunc){
	x->re+=(pow(couplings[KAPPA_A],2)+pow(couplings[KAPPA_S],2))*
	  elfkappa(left[0]);
      }
      /* Half the contribution of any gluons */
      if(ngl>=1)
	x->re+=0.25*couplings[MM]/GE(left[1]);
    }
  }

    /* qqgg gauge interactions */

    if(npl==npr && ngl>=1){
      for(j=2; j<npl; j++)if(left[j]!=right[j])break;
      if(j==npl){
	if(left[0]==right[0] && H_EQUAL(left[1],right[1]))
	  x->re+=couplings[GGN]*elfquarkglue(left[0],left[1]);
	else if(H_EQUAL(left[0],right[0]) && H_EQUAL(left[1],right[1]))
	  x->re-=couplings[GGN]*(GE(left[1])+GE(right[1]))/
	    (8.0*PI*sqrt(GE(left[1])*GE(right[1]))*
	     pow(QE(left[0])-QE(right[0]),2));
	
      }
    }
    else if(npl+2==npr){
      for(j=1; j<npl; j++)if(left[j]!=right[j+2])break;
      if(j==npl){
	if(H_EQUAL(left[0],right[0]) && H_MINUS(right[1],right[2])){
	   temp=-couplings[GGN]*(GE(right[1])-GE(right[2]))/
	   (8*PI*sqrt(GE(right[1])*GE(right[2]))*
	    pow(QE(left[0])-QE(right[0]),2));
	   x->re+=temp;
	}
      }
    }
    else if(npr+2==npl){
      for(j=1; j<npr; j++)if(left[j+2]!=right[j])break;
      if(j==npr){
	if(H_EQUAL(left[0],right[0]) && H_MINUS(left[1],left[2]))
	  x->re-=couplings[GGN]*(GE(left[1])-GE(left[2]))/
	    (8*PI*sqrt(GE(left[1])*GE(left[2]))*
		       pow(QE(left[0])-QE(right[0]),2));
      }
    }
    
   /* qqg Yukawa interactions */
    
  if(npl+1==npr){
    for(j=1; j<npl; j++)if(left[j]!=right[j+1])break;
    if(j==npl){
      if(H_EQUAL(left[0],right[0])){
	x->re-=couplings[KAPPA_S]*couplings[MU_F]*(QE(left[0])+QE(right[0]))/
	  (4*sqrt(PI*GE(right[1]))*QE(left[0])*QE(right[0]));
      }
      temp=couplings[KAPPA_A]*couplings[MU_F]*sqrt(GE(right[1]))/
	(4*sqrt(PI)*QE(right[0])*QE(left[0]));
      ZTIMES(x,&zeta(left[0],right[1],right[0]),temp);
    }
  } 
  else if(npr+1==npl){
    for(j=1; j<npr; j++)if(left[j+1]!=right[j])break;
    if(j==npr){
      if(H_EQUAL(left[0],right[0])){
	x->re-=couplings[KAPPA_S]*couplings[MU_F]*
	  (QE(left[0])+QE(right[0]))/
	  (4*sqrt(PI*GE(left[1]))*QE(left[0])*QE(right[0]));
      }
      temp=couplings[KAPPA_A]*couplings[MU_F]*
	sqrt(GE(left[1]))/(4*sqrt(PI)*QE(right[0])*QE(left[0]));
      ZTIMES(x,&zeta(left[0],left[1],right[0]),temp);
    }
  }
  
  /* qqgg Yukawa interactions */

#if YUKAWA_FOUR_POINT  /* flag to include meson 4-point Yukawa interactions */ 
  if(npl==npr && ngl>=1){
    for(j=2; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      if(H_EQUAL(left[0],right[0]))
	x->re+=pow(couplings[KAPPA_S],2)/
	  (8*PI*sqrt(GE(left[1])*GE(right[1]))*(QE(right[0])+GE(right[1])));
      temp=couplings[KAPPA_A]*couplings[KAPPA_S]/
	(8*PI*sqrt(GE(left[1])*GE(right[1]))*(QE(right[0])+GE(right[1])));
      ZTIMES(x,&zeta(left[0],left[1],right[0]),temp);
      ZTIMES(x,&zeta(left[0],right[1],right[0]),temp);
      temp=pow(couplings[KAPPA_A],2)/
	(8*PI*sqrt(GE(left[1])*GE(right[1]))*(QE(right[0])+GE(right[1])));
      ZTIMES(x,&zetazeta(left[0],left[1],right[1],right[0]),temp);
    }
  }
  else if(npl+2==npr){
    for(j=1; j<npl; j++)if(left[j]!=right[j+2])break;
    if(j==npl){
      if(H_EQUAL(left[0],right[0]))
	x->re+=pow(couplings[KAPPA_S],2)/
	  (8*PI*sqrt(GE(right[1])*GE(right[2]))*(QE(right[0])+GE(right[1])));
      temp=couplings[KAPPA_A]*couplings[KAPPA_S]/
	(8*PI*sqrt(GE(right[1])*GE(right[2]))*(QE(right[0])+GE(right[1])));
      ZTIMES(x,&zeta(left[0],right[1],right[0]),temp);
      ZTIMES(x,&zeta(left[0],right[2],right[0]),-temp);
      temp=-pow(couplings[KAPPA_A],2)/
	(8*PI*sqrt(GE(right[1])*GE(right[2]))*(QE(right[0])+GE(right[1])));
      ZTIMES(x,&zetazeta(left[0],right[2],right[1],right[0]),temp);
    }
  }
  else if(npr+2==npl){
    for(j=1; j<npr; j++)if(left[j+2]!=right[j])break;
    if(j==npr){
      if(H_EQUAL(left[0],right[0]))
	x->re+=pow(couplings[KAPPA_S],2)/
	  (8*PI*sqrt(GE(left[1])*GE(left[2]))*
	   (QE(left[0])+GE(left[1])));
      temp=couplings[KAPPA_A]*couplings[KAPPA_S]/
	(8*PI*sqrt(GE(left[1])*GE(left[2]))*(GE(left[2])-QE(right[0])));
      ZTIMES(x,&zeta(left[0],left[2],right[0]),temp);
      ZTIMES(x,&zeta(left[0],left[1],right[0]),-temp);
      temp=pow(couplings[KAPPA_A],2)/
	(8*PI*sqrt(GE(left[1])*GE(left[2]))*(GE(left[2])-QE(right[0])));
      ZTIMES(x,&zetazeta(left[0],left[1],left[2],right[0]),temp);
    }
  }
#endif

  return;
}

/*
   quark-gluon interactions with nonzero momentum included.
   quark-antiquark interactions are not included
   x is not initialized
   ngl is number of gluons in left state.
*/

void interact_quarkc(complex_element *x, partype *left, int npl, 
		     partype *right, int npr, element *couplings, int ngl){
  int j;
  element temp;
  complex_element ztemp,z;
  

  /* Mass terms and Yukawa self-intertias 
     we don't have any phase factors here */

  if(npl==npr){
    for(j=0; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      temp=pow(couplings[MU_F],2);
      if(npl==working.nptrunc)
	temp+=couplings[MMF_0];
      if(npl==working.nptrunc-1)
	temp+=couplings[MMF_1];
      if(npl<=working.nptrunc-2)
	temp+=couplings[MMF_2];
       x->re+=temp/(2.0*QE(left[0]));
      if(npl<working.nptrunc){
	x->re+=(pow(couplings[KAPPA_A],2)+pow(couplings[KAPPA_S],2))*
	  elfkappa(left[0]);
      }
      /* Half the contribution of any gluons */
      if(ngl>=1)
	x->re+=0.25*couplings[MM]/GE(left[1]);
    }
  }

    /* qqgg gauge interactions */

  if(npl==npr && ngl>=1){
    for(j=2; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      quark_phase(&ztemp,left,2,right,2);	
      if(left[0]==right[0] && H_EQUAL(left[1],right[1]))
	x->re+=couplings[GGN]*elfquarkglue(left[0],left[1]);
      else if(H_EQUAL(left[0],right[0]) && H_EQUAL(left[1],right[1])){
	temp=-couplings[GGN]*(GE(left[1])+GE(right[1]))/
	  (8.0*PI*sqrt(GE(left[1])*GE(right[1]))*
	   pow(QE(left[0])-QE(right[0]),2));
	ZTIMES(x,&ztemp,temp);
      }
    }
  }
  else if(npl+2==npr){
    for(j=1; j<npl; j++)if(left[j]!=right[j+2])break;
    if(j==npl){
      quark_phase(&ztemp,left,1,right,3);
      if(H_EQUAL(left[0],right[0]) && H_MINUS(right[1],right[2])){
	temp=-couplings[GGN]*(GE(right[1])-GE(right[2]))/
	  (8*PI*sqrt(GE(right[1])*GE(right[2]))*
	   pow(QE(left[0])-QE(right[0]),2));
	ZTIMES(x,&ztemp,temp);
      }
    }
  }
  else if(npr+2==npl){
    for(j=1; j<npr; j++)if(left[j+2]!=right[j])break;
    if(j==npr){
      quark_phase(&ztemp,left,3,right,1);
      if(H_EQUAL(left[0],right[0]) && H_MINUS(left[1],left[2])){
	temp=-couplings[GGN]*(GE(left[1])-GE(left[2]))/
	  (8*PI*sqrt(GE(left[1])*GE(left[2]))*
		     pow(QE(left[0])-QE(right[0]),2));
	ZTIMES(x,&ztemp,temp);
      }
    }
  }
    
    /* qqg Yukawa interactions */

  if(npl+1==npr){
    for(j=1; j<npl; j++)if(left[j]!=right[j+1])break;
    if(j==npl){
      quark_phase(&z,left,1,right,2);
      if(H_EQUAL(left[0],right[0])){
	temp=-couplings[KAPPA_S]*couplings[MU_F]*(QE(left[0])+QE(right[0]))/
	  (4*sqrt(PI*GE(right[1]))*QE(left[0])*QE(right[0]));
	ZTIMES(x,&z,temp);
      }
      temp=couplings[KAPPA_A]*couplings[MU_F]*sqrt(GE(right[1]))/
	(4*sqrt(PI)*QE(right[0])*QE(left[0]));
      ZZPROD(&ztemp,&zeta(left[0],right[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
    }
  } 
  else if(npr+1==npl){
    for(j=1; j<npr; j++)if(left[j+1]!=right[j])break;
    if(j==npr){
      quark_phase(&z,left,2,right,1);
      if(H_EQUAL(left[0],right[0])){
	temp=-couplings[KAPPA_S]*couplings[MU_F]*
	  (QE(left[0])+QE(right[0]))/
	  (4*sqrt(PI*GE(left[1]))*QE(left[0])*QE(right[0]));
	ZTIMES(x,&z,temp);
      }
      temp=couplings[KAPPA_A]*couplings[MU_F]*
	sqrt(GE(left[1]))/(4*sqrt(PI)*QE(right[0])*QE(left[0]));
      ZZPROD(&ztemp,&zeta(left[0],left[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
    }
  }
  
  /* qqgg Yukawa interactions */

#if YUKAWA_FOUR_POINT  /* flag to include meson 4-point Yukawa interactions */ 
  if(npl==npr && ngl>=1){
    for(j=2; j<npl; j++)if(left[j]!=right[j])break;
    if(j==npl){
      quark_phase(&z,left,2,right,2);
      if(H_EQUAL(left[0],right[0])){
	temp=pow(couplings[KAPPA_S],2)/
	  (8*PI*sqrt(GE(left[1])*GE(right[1]))*(QE(right[0])+GE(right[1])));
	ZTIMES(x,&z,temp);
      }
      temp=couplings[KAPPA_A]*couplings[KAPPA_S]/
	(8*PI*sqrt(GE(left[1])*GE(right[1]))*(QE(right[0])+GE(right[1])));
      ZZPROD(&ztemp,&zeta(left[0],left[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
      ZZPROD(&ztemp,&zeta(left[0],right[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
      temp=pow(couplings[KAPPA_A],2)/
	(8*PI*sqrt(GE(left[1])*GE(right[1]))*(QE(right[0])+GE(right[1])));
      ZZPROD(&ztemp,&zetazeta(left[0],left[1],right[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
    }
  }
  else if(npl+2==npr){
    for(j=1; j<npl; j++)if(left[j]!=right[j+2])break;
    if(j==npl){
      quark_phase(&z,left,1,right,3);
      if(H_EQUAL(left[0],right[0])){
	temp=pow(couplings[KAPPA_S],2)/
	  (8*PI*sqrt(GE(right[1])*GE(right[2]))*(QE(right[0])+GE(right[1])));
	ZTIMES(x,&z,temp);
      }
      temp=couplings[KAPPA_A]*couplings[KAPPA_S]/
	(8*PI*sqrt(GE(right[1])*GE(right[2]))*(QE(right[0])+GE(right[1])));
      ZZPROD(&ztemp,&zeta(left[0],right[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
      ZZPROD(&ztemp,&zeta(left[0],right[2],right[0]),&z);
      ZTIMES(x,&ztemp,-temp);
      temp=-pow(couplings[KAPPA_A],2)/
	(8*PI*sqrt(GE(right[1])*GE(right[2]))*(QE(right[0])+GE(right[1])));
      ZZPROD(&ztemp,&zetazeta(left[0],right[2],right[1],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
    }
  }
  else if(npr+2==npl){
    for(j=1; j<npr; j++)if(left[j+2]!=right[j])break;
    if(j==npr){
      quark_phase(&z,left,3,right,1);
      if(H_EQUAL(left[0],right[0])){
	temp=pow(couplings[KAPPA_S],2)/
	  (8*PI*sqrt(GE(left[1])*GE(left[2]))*(QE(left[0])+GE(left[1])));
	ZTIMES(x,&z,temp);
      }
      temp=couplings[KAPPA_A]*couplings[KAPPA_S]/
	(8*PI*sqrt(GE(left[1])*GE(left[2]))*(GE(left[2])-QE(right[0])));
      ZZPROD(&ztemp,&zeta(left[0],left[2],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
      ZZPROD(&ztemp,&zeta(left[0],left[1],right[0]),&z);
      ZTIMES(x,&ztemp,-temp);
      temp=pow(couplings[KAPPA_A],2)/
	(8*PI*sqrt(GE(left[1])*GE(left[2]))*(GE(left[2])-QE(right[0])));
      ZZPROD(&ztemp,&zetazeta(left[0],left[1],left[2],right[0]),&z);
      ZTIMES(x,&ztemp,temp);
    }
  }
#endif
  
  return;
}