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

/*
    This file contains subroutines to calculate improved 
    matrix elements for DLCQ.  

    This calculates P^- (previously, it gave 2 P^-, this was changed in 
    Feb. 2002, as noted in comments)

    This follows the basic method
    of Appendix C of hep-ph/9704408 with the following exceptions:

    1. The method is now numerically stable with values good
    to machine precision.  This is because a basis of sines and
    cosines is used instead of the Polynomials.

    3. For particle number changing interactions, the wavefunction
    is assumed to go to a constant as any two momenta go to
    zero (this can be proven from the bound state equations).

    4. There is no limit to K other than K = KKMAX.

    The method for calculating the matrix elements of the instantaneous 
    interaction rely  on a combination of numerical Gauss-Jacobi
    integration (The number of necessary grid points was determined
    before hand in Mathematica) and analytic integration.

    External subroutines to provide the weighting functions,
    the gamma function, and Bessel functions are used.        
                                                                     
   In the parallelized version, all order K^2 computations
   are parallelized.  This includes the numerical integrations
   in instant22() the matrix multiplications in mapcouplings();

   These are all performed with the communications handle
   IMPROVED_MPI_COMM;

   For large K, the arrays holding the matrix elements can
   be very large, increasing as K^3.  The use of order_k(...) 
   minimizes the storage needed for the 2->2 particle interactions. 
*/

#ifdef USE_MPI
#include "usempi.h"

#if 0 /* GROUP_SIZE != 0*/   /* case with master and slaves */
#define IMPROVED_MPI_COMM mpiallslaves
#else  /* case where everyone does the work */
#define IMPROVED_MPI_COMM MPI_COMM_WORLD
#endif

#endif

/* 
   mapcouplings(...) Calculate beta and calculate matrix elements
   in the arrays fun22 and fun13.
   ratio=m^2=mu^2 a^{D-2}/g^2 N and kt is twice the 
   maximum value of K to calculate matrix elements.
   The root finding for beta is done 
   internally since efficiency is not very important here.
   

   Here, I have rescaled g^2 -> g^2/2.0 so that the definition
   of g is now standard.   
   
   If the value of beta matches the previous value and
   global_kt is the largest value of kt for any spectrum,
   as determined by initspectrum(...).
*/

unsigned int global_kt=0;  /* initial value for global_kt */
/* the following variables are used to determine
   if new matrix elements need to be calculated */
static element last_beta=-1.0;    /* previous value of beta */
static unsigned int last_kt=0; /* previous value of kt   */
size_t fun22_memory=0; /* bytes of memory used fun22 currently used */
size_t fun13_memory=0; /* bytes of memory used fun13 currently used */

/*  This file uses a lot of functions which are only
    needed locally */

/*  This is the number terms in fourelements and two elements */
const int four=4,two=2;

#define TOP (3*KKMAX+20) /*This is the maximum number of states */
                         /*in the cosine and sine bases */ 
                         /*This definition should match the definition
                           of top below. */

#if IMPROVED==3
static struct {element x; element b; element tol; integer which; 
               element k1; element k2; element u; element kk;} f;
static element bj(element, element);
static void inner22(element *,element *,int,element,element);
static void map22(element *,element *,int,int,element,element);
static void mass22(element *,element *,int,element,element);
static void instant22(element *,element *,int,element,element);
static void annihilate22(element *,element *,int,element,element);
static void lambda22(element *,element *,int,element,element);
#endif

element mapcouplings(element *couplings, unsigned int this_kt){
  int i,j,k;
  element beta=0.0,pi,ratio=couplings[MM]/couplings[GGN];
  element lower=0.0, upper=0.5, tol=10.0*DBL_EPSILON;
  clock_t t1,t2;
  size_t size_temp;
  integer kkdown; /* sum of two incoming momenta, *each* rounded down */
#if IMPROVED==1
  element k1,k2,k3,k4,k123,kl;
  int l;
  element ksum;
#elif IMPROVED==3
  element temp;
  fun13_forms *sow2;
  int l,m;
  integer kk2,top,error,itype=1,lwork;
  element *vvi[2],*vvl[2],*ham[2],*tt[2],*tt2;
  element offset=-0.5,pre;
  char up='U',vects='V';
#else
  element *rowin;
#endif
  fun22_forms *row;
  fun13_forms *sow;
#ifdef USE_MPI
  /* initialize communicator with just slaves. */
  int nslaves,myrank;
#endif
  pi= atan(1.) * 4.;

  /* Complain if this_kt is larger than global_kt */

  if(this_kt>global_kt){
    fprintf(stderr,__FILE__":  present value of kt=%i is larger than "
	    "global_kt=%i\n   Use initspectrum() to properly initialize "
	    "global_kt.\n",this_kt,global_kt);
    goto error;
  }

  /* Complain if global_kt is too large */

  if(global_kt>2*KKMAX){
     fprintf(stderr,__FILE__":  Value of kt=%i in mapcouplings is too large,"
	    " exiting\n",global_kt);
    goto error;
  }
  
  /* Complain if improved and periodic */

  if(GLUE_PERIODIC && IMPROVED==3){
    fprintf(stderr,__FILE__":  Improved matrix elements for periodic\n"
	    "     boundary conditions does not work yet.\n");
    goto error;
  }

/*  Calculate beta and prepare coupings. */
  if(ratio<tol)
    beta=0.0;
  else {
    while(upper-lower>tol){
      beta=0.5*(upper+lower);
      if(beta * tan(pi * beta)<ratio) /* factor of 2.0 removed */
	lower=beta;
      else
	upper=beta;
    }
  }

  /* see if anything new needs to be calculated,
     if not, just return                            */

  if(fabs(beta-last_beta)<tol && global_kt<=last_kt){
    return beta;
  }

  last_kt=global_kt;
  last_beta=beta;
  t1=t2=clock();
#ifdef USE_MPI
  MPI_Comm_size(IMPROVED_MPI_COMM,&nslaves);
  MPI_Comm_rank(IMPROVED_MPI_COMM,&myrank);
#endif
  
  /* calculate 2->2 particle interactions */
#if 0 /*debugging print */ 
  printf("Starting 2->2 max K=%i/2 beta=%f\n",global_kt,beta);
#endif
  
#if IMPROVED==3 /* internally constructed improved matrix elements */
  /* top is the number of wavefunctions to include in the basis;
     it should be roughly proportional to K.
     The time needed to calculate the improved matrix elements
     is roughly order top^4.                             */
  /* Timing measured on rum or whiskey  */
  /* top=4*K, K=20 takes 140 seconds.  matrix elements 2.5 digits. */
  /* top=3*K, K=20 takes 54 seconds.  matrix elements 2.5 digits.  */
  /* top=2*K, K=20 takes 15 seconds. matrix elements 2 digits.  */
  top=(GLUE_PERIODIC?3*global_kt+40:(3*global_kt+40)/2);
  if(top>TOP){
    fprintf(stderr,__FILE__":  Value of top=%li in mapcouplings is too large,"
	    " exiting\n",top);
    goto error;
  }
  for(i=0; i<2; i++){
    if(NULL==(vvi[i]=malloc(top*top*sizeof(element))))goto memerror;
    if(NULL==(vvl[i]=malloc(top*top*sizeof(element))))goto memerror;
    if(NULL==(ham[i]=malloc(top*top*sizeof(element))))goto memerror;
    if(NULL==(tt[i]=malloc(top*(global_kt+3)/2*
			   sizeof(element))))goto memerror;
  }
  if(NULL==(tt2=malloc(top*(global_kt+3)/2*
		       sizeof(element))))goto memerror;
  t2=clock();
  instant22(vvi[0],vvi[1],top,beta,offset);
  t2=clock()-t2;
  lambda22(vvl[0],vvl[1],top,beta,offset);
  if(beta>tol){
    mass22(ham[0],ham[1],top,beta,offset);
    pre=beta*tan(pi*beta);
    for(k=0; k<2; k++)
      for(i=0; i<top; i++)
	for(j=0; j<=i; j++)
	  ham[k][i*top+j]=2.0*(vvi[k][i*top+j]+
			       pre*ham[k][i*top+j])-vvl[k][i*top+j];
  } else {
    for(k=0; k<2; k++)
      for(i=0; i<top; i++)
	for(j=0; j<=i; j++)
	  ham[k][i*top+j]=2.0*vvi[k][i*top+j]-vvl[k][i*top+j];
  }
  inner22(vvl[0],vvl[1],top,beta,offset);
#if 0 /* print out hamiltonian and inner product matrices */
  for(k=0; k<2; k++)
    for(i=0; i<top; i++)
      for(j=0; j<=i; j++)
	printf("  %i %i %i %.16e %.16e \n",
	       k,i,j,ham[k][i*top+j],vvl[k][i*top+j]);
#endif
  /* Solve the generalized eigenvalue problem */
  /* The eigenvectors are outputted in ham */
  /* Can only access one term of tt[0]  */
  lwork=top*((global_kt+3)/2);
  for(k=0; k<2; k++){
    DSYGV(&itype,&vects,&up,&top,ham[k],&top,vvl[k],&top,
	  tt[0],tt[1],&lwork,&error,1,1);
    if(error){
      fprintf(stderr,__FILE__":  Error=%li in DSYGV in mapcouplings\n",
	      error);
      fprintf(stderr,__FILE__":  k=%i,itype=%li,vects=%c,up=%c,top=%li\n",
	      k,itype,vects,up,top);
      goto error;
    }
#if 0 /* print out eigenvalues and eigenvectors. */
    printf("eigenvalues and eigenvectors, k=%i\n",k);
    for(i=0; i<top; i++){
      printf(" %.16e:  ",tt[0][i]);
      for(j=0; j<top; j++)
	printf(" %.6e",ham[k][i*top+j]);
      printf("\n");
    }
#endif
  }
  
  
  /*  Loop through different sectors of K and make improved
      matrix elements */
  /* kkdown is the sum of the momenta, *each* rounded down */
  for(kkdown=0; kkdown<=global_kt/2-(GLUE_PERIODIC?0:1); kkdown++){
    /* kk is the sum of the two incoming momenta */
    integer kk=(GLUE_PERIODIC?kkdown:kkdown+1);
    /* don't try zero mode sector */
    if(kk==0)continue;
    /* j is one incoming momentum (rounded down) */
    for(j=0; j<=kkdown; j++){
      if(fun22[kkdown][j]==NULL){
	size_temp=((2*j<=kkdown?j:kkdown-j)+1)*sizeof(fun22_forms);
        fun22_memory+=size_temp;
	if(NULL==(fun22[kkdown][j]=malloc(size_temp)))
	  goto memerror;
      }
    }
    
#ifdef USE_MPI
    if(kkdown%nslaves==myrank)
#endif
      {
	map22(tt[0],tt[1],top,kk,beta,offset);
#if 0 /* debugging print of mapping matrix*/
	printf("mapping matrix, kk=%li\n",kk);
	for(k=0; k<2; k++)
	  for(j=0; j<kk; j++){
	    for(i=0; i<top; i++)
	      printf(" %.16e\n",tt[k][i*kk+j]);
	    printf("\n");
	  }
#endif
	/* the second half of tt2[] is reserved for the signs */
	lwork=top*top-kk;
	for(k=0; k<2; k++){
	  kk2=(kk+1-k)/2;
	  for(j=0; j<kk; j++)
	    for(l=0; l<kk2; l++){
	      tt2[l*kk+j]=0.0;
	      for(i=0; i<top; i++){
		tt2[l*kk+j]+=tt[k][i*kk+j]*ham[k][l*top+i];
	      }
	    }
	  /* Perform QR factorization */
	  /* Unfortunately, R-often does not have positive diagonals,
	     thus, we need to switch signs.  
	     This is essential for good performance */
	  DGEQRF(&kk,&kk2,tt2,&kk,vvl[0],vvl[0]+kk,&lwork,&error);
	  if(error){
	    fprintf(stderr,__FILE__":  DGEQRF error %li in mapcouplings\n",
		    error);
	    goto error;
	  }
	  for(l=0; l<kk2; l++)
	    if(tt2[l*kk+l]>0.0)
	      vvl[1][l]=1.0; 
	    else 
	      vvl[1][l]=-1.0;
	  DORGQR(&kk,&kk2,&kk2,tt2,&kk,vvl[0],vvl[0]+kk,&lwork,&error);
	  if(error){
	    fprintf(stderr,__FILE__":  DORGQR error %li in mapcouplings\n",
		    error);
	    goto error;
	  }
	  /* Perform multiplication to get to cosine-sine basis */
	  for(j=0; j<kk; j++)
	    for(i=0; i<top; i++){
	      tt[k][i*kk+j]=0.0;
	      for(l=0; l<kk2; l++)
		tt[k][i*kk+j] += vvl[1][l]*tt2[l*kk+j]*ham[k][l*top+i];
	    }
	}
#if 0 /* print out transformation matrix */
	printf("transform matrix, kk=%li\n",kk);
	printf("lwork=%li, top=%li\n",lwork,top);
	for(k=0; k<2; k++)
	  for(j=0; j<kk; j++){
	    for(i=0; i<top; i++)
	      printf(" %.16e\n",tt[k][i*kk+j]);
	    printf("\n");
	  }
#endif
	
	/* Loop through different interaction terms */
	for(m=0; m<four; m++){
	  switch(m){
	  case 0:
	    mass22(vvl[0],vvl[1],top,beta,offset);
	    break;
	  case 1:
	    for(l=0; l<2; l++)
	      memcpy(vvl[l],vvi[l],top*top*sizeof(doublereal));
	    break;
	  case 2:
	    annihilate22(vvl[0],vvl[1],top,beta,offset);
	    break;
	  case 3:
	    lambda22(vvl[0],vvl[1],top,beta,offset);
	    break;
	  }
	  /* Perform matrix multiplications to get back to DLCQ */
	  for(k=0; k<2; k++){
	    for(j=0; j<=kkdown; j++)
	      for(i=0; i<top; i++){
		tt2[i*(kkdown+1)+j]=0.0;
		for(l=0; l<top; l++)
		  tt2[i*(kkdown+1)+j]+=
		    tt[k][l*(kkdown+1)+j]*vvl[k][l*top+i];
	      }
	    for(j=0; j<=kkdown; j++)
	      for(i=0; i<=kkdown; i++){
		vvl[k][i*(kkdown+1)+j]=0.0;
		for(l=0; l<top; l++)
		  vvl[k][i*(kkdown+1)+j]+=
		    tt[k][l*(kkdown+1)+j]*tt2[l*(kkdown+1)+i];
	      }
	  }
	  /* put results into fun22 */
	  for(j=0; j<=kkdown; j++)
	    for(k=0; k<=j && k<=kkdown-j; k++){
	      temp=(vvl[0][j*(kkdown+1)+k]+
		    vvl[1][j*(kkdown+1)+k])/(element) kk;
	      row=&fun22[kkdown][j][k];
	      switch(m){
	      case 0:
		row->mass=temp;
		break;
	      case 1:
		row->ggn=temp;
		break;
	      case 2:
		row->annihilate=temp;
		break;
	      case 3:
		row->lambda=temp;
		break;
	      }
	    }
	}
      }
  }
  
  /* gather and send out results */
  
#ifdef USE_MPI
  /* kkdown is the sum of the momenta, *each* rounded down */
  for(kkdown=0; kkdown<=global_kt/2-(GLUE_PERIODIC?0:1); kkdown++)
    for(j=0; j<=kkdown; j++)
      MPI_Bcast(fun22[kkdown][j],((2*j<=kkdown?j:kkdown-j)+1)*
		sizeof(fun22_forms)/sizeof(element),
		MPI_DOUBLE,
		kkdown%nslaves,IMPROVED_MPI_COMM);
#endif
  
  /* free in reverse order of allocation */
  
  free(tt2);
  for(i=1; i>=0; i--){
    free(tt[i]); free(ham[i]); free(vvl[i]); free(vvi[i]);
  }
  
#else  /* internally constructed DLCQ */
  /* kkdown is the sum of the momenta, *each* rounded down */
  for(kkdown=0; kkdown<=global_kt/2-(GLUE_PERIODIC?0:1); kkdown++){
    /* j is one incoming momentum (rounded down) */
    for(j=0; j<=kkdown; j++){
      if(fun22[kkdown][j]==NULL){
	size_temp=((2*j<=kkdown?j:kkdown-j)+1)*sizeof(fun22_forms);
	fun22_memory+=size_temp;
	if(NULL==(fun22[kkdown][j]=malloc(size_temp)))
	  goto memerror;
      }
      /* k is one outgoing momentum (rounded down) */
      for(k=0; k<=j && k<=kkdown-j; k++){
	row=&fun22[kkdown][j][k];
	ksum=(element) (GLUE_PERIODIC?kkdown:kkdown+1);
	k1=GE(j); k2=ksum-k1;
	k3=GE(k); k4=ksum-k3;
	if(j==k)
	  row->mass=0.25*(1.0/k1+1.0/k2);   /* changed 0.5 -> 0.25 */
	else 
	  row->mass=0.0;
	if(j==k){
	  row->ggn=1.0/(16.0*sqrt(k1*k2)); /* changed 8->16 */
	  for(l=0; l<=kkdown; l++){
	    /* no zero mode */
	    if(GLUE_PERIODIC && (l==0 || l==kkdown))continue; 
	    if(l==j)continue; /* no pole */
	    kl=GE(l);
	    row->ggn+=(k1+kl)*(k2+ksum-kl)/
	      (16*pi*sqrt(k1*k2*kl*(ksum-kl))*pow(k1-kl,2)); /* 8->16 */
	  }
	} else
	  row->ggn=-(k1+k3)*(k2+k4)/
	    (16.0*pi*sqrt(k1*k2*k3*k4)*pow(k1-k3,2));         /* 8->16 */  
	row->annihilate=-(ksum-2.0*k1)*(ksum-2.0*k3)/
	  (16.0*pi*sqrt(k1*k2*k3*k4)*ksum*ksum);              /* 8->16 */
	row->lambda=1.0/(8.0*pi*sqrt(k1*k2*k3*k4)); /* 4.0 -> 8.0 */
      }
    }
  }
#endif /* value of improved */
  
  /**** calculate 1->3 particle interactions *****/
#if 0 /*debugging print */ 
  printf("Starting 1->3\n");
#endif
  /* i is the momentum (rounded down) of the single incoming particle */
  for(i=(GLUE_PERIODIC?0:1); i<=(global_kt-(GLUE_PERIODIC?0:1))/2; i++){
    /* j is the momentum (rounded down) of an outgoing particle */
    for(j=0; j<=(GLUE_PERIODIC?i:i-1); j++){ 
      
      if(fun13[i][j]==NULL){
	size_temp=(i-j+1)*sizeof(fun13_forms);
	fun13_memory+=size_temp;
	if(NULL==(fun13[i][j]=malloc(size_temp)))
	  goto memerror;
      }
      
      
#if IMPROVED==3  /* construct improved matrix elements */
#ifdef USE_MPI
      if(i%nslaves==myrank)
#endif
	/* k is momentum (rounded down) of another outgoing particle */
	for(k=0; k+j<=(GLUE_PERIODIC?i:i-1) && k<=j; k++){
	  /* don't do zero modes */
	  if(GLUE_PERIODIC && (k==0 || j==0 || i==k+j))continue; 
	  sow=&fun13[i][j][k];
	  temp=1.0/(8.0*pi*sqrt(GE(i)*
				inner13(1,i,j,k,beta,&error)));  /* 4->8 */
	  sow->ggn=0.5*temp*inner13(3,i,j,k,beta,&error);
	  sow->lambda=temp*inner13(2,i,j,k,beta,&error);
	  if(error){
	    fprintf(stderr,__FILE__":  inner13 error=%li in mapcouplings\n",
		    error);
	    goto error;
	  }
	  if(k!=j){	      
	    sow2=&fun13[i][k][j];
	    sow2->ggn =- sow->ggn;
	    sow2->lambda = sow->lambda;
	  }
	}
      
#else      /* construct DLCQ matrix elements */
      for(k=0; k+j<=(GLUE_PERIODIC?i:i-1); k++){
	sow=&fun13[i][j][k];
	k123=GE(i);
	k1=GE(j);
	k2=GE(k);
	k3=k123-k1-k2;
	sow->ggn=(k1-k2)*(k123+k3)/
	  (16.0*pi*sqrt(k1*k2*k3*k123)*pow(k1+k2,2));  /* 8->16 */
	sow->lambda=1.0/(8.0*pi*sqrt(k1*k2*k3*k123));       /* 4->8 */
      }
#endif
    }
  }
  
  /* gather and send out results */
  
#if defined(USE_MPI) && IMPROVED==3
  for(i=(GLUE_PERIODIC?0:1); i<=(global_kt-(GLUE_PERIODIC?0:1))/2; i++)
    /* j is the momentum (rounded down) of an outgoing particle */
    for(j=0; j<=(GLUE_PERIODIC?i:i-1); j++)
      MPI_Bcast(fun13[i][j],(i-j+1)*sizeof(fun13_forms)/sizeof(element),
		MPI_DOUBLE,
		i%nslaves,IMPROVED_MPI_COMM);
#endif
  
#if 0 /*********debugging print *****************************/
  printf("finished in a total of %f seconds (instant in %f).\n    ",
	 (float)(clock()-t1)/(float) CLOCKS_PER_SEC,
	 (float) t2/(float) CLOCKS_PER_SEC); 
  /* kkdown is the sum of the momenta, *each* rounded down */
  for(kkdown=0; kkdown<5 && kkdown<=global_kt/2-(GLUE_PERIODIC?0:1); 
      kkdown++){
    int kk=(GLUE_PERIODIC?kkdown:kkdown+1);
    double x=1.0;
    printf("starting Ktot=%i%s times %f\n",kk,(GLUE_PERIODIC?"":"/2"),x);
    for(j=0; j<=kkdown; j++)
      for(k=0; k<=j && k<=kkdown-j; k++){
	row=&fun22[kkdown][j][k];
	printf("%.12f %.12f %.12f %.12f\n",x*row->mass,x*row->ggn,
	       x*row->annihilate,x*row->lambda);
      }
  } 
  /* i is the momentum (rounded down) of the single incoming particle */
  for(i=(GLUE_PERIODIC?0:1); i<5 && i<=(global_kt-(GLUE_PERIODIC?0:1))/2; i++){
    double x=2.0;
    printf("starting Kin=%i%s times %f\n",GLUE_PERIODIC?i:2*i+1,
	     GLUE_PERIODIC?"":"/2",x);
    /* j is the momentum (rounded down) of an outgoing particle */
    for(j=0; j<=(GLUE_PERIODIC?i:i-1); j++)
      for(k=0; k+j<=(GLUE_PERIODIC?i:i-1); k++){
	sow=&fun13[i][j][k];
	printf("%.12f %.12f\n",x*sow->ggn,x*sow->lambda);
      }
  }
#endif /* end debugging prints */
  
#ifdef USE_MPI /* all slaves must finish before going further. */
  MPI_Barrier(IMPROVED_MPI_COMM);
#endif
  
  return beta; /* normal return */

memerror:
   fprintf(stderr,__FILE__":  malloc returning null in mapcouplings\n");
error:
   fprintf(stderr,__FILE__":  Error in mapcouplings\n");
#ifdef USE_MPI
   fprintf(stderr,__FILE__":  Slave %i sending abort signal\n",myrank);
  MPI_Abort(IMPROVED_MPI_COMM,1);
#endif
  return -1.0;
}

/******************************************************************/
/******************************************************************/
/******************************************************************/


#if IMPROVED==3 /* Functions for internal improved matrix element */

/* The function J_n(x/2)/x^n  */

static element bj(element n, element x){
  doublereal n1,x1;
  if(n<0.0)
    return 4.0*(n+1.0)*bj(n+1.0,x)-x*x*bj(n+2.0,x);
  else if(fabs(x) < ELEMENT_EPSILON){
    n1=n+1;
    return 1.0/(MYGAMMA(&n1)*pow(4.0,n));
  }
  else {
    n1=n; x1=0.5*fabs(x);
    return MYBESSELJ(&n1,&x1)/pow(fabs(x),n);
  }
}


/*  Inner produces top by top matrices for cosine and sine bases
    One triangle of the matrix is filled. */

static void inner22(element *cc, element *ss, int top, 
		    element beta, element offset){
  int i,j;
  element pre,pi,m,n;
  doublereal b2=2.0*beta+1.0;
  pi=4.0*atan(1.0);
  pre=0.5*sqrt(pi)*MYGAMMA(&b2);
  for(i=0; i<top; i++)
    for(j=0; j<=i; j++){
      m=2.0*(element) i+offset; n=2.0*(element) j+offset;
      cc[i*top+j]=pre*(bj(2*beta+0.5,pi*(m-n))+
			  bj(2*beta+0.5,pi*(m+n)));
      m += 1.0; n += 1.0;
      ss[i*top+j]=pre*(bj(2*beta+0.5,pi*(m-n))-
			  bj(2*beta+0.5,pi*(m+n)));
    }
}

/* translate inner product basis to DLCQ basis
   coscos is top by kk
   sinsin to top by kk              */

static void map22(element *cc, element *ss, int top, 
		  int kk, element beta, element offset){
  int i,j;
  element c1,pi,x,m;
  pi=4.0*atan(1.0);
  for(i=0; i<top; i++)
    for(j=0; j<kk; j++){
      x=(0.5+(element) j)/((element) kk);
      c1=pow(x*(1.0-x),beta);
      m=2.0*(element)i+offset;
      cc[i*kk+j]=c1*cos(m*pi*(x-0.5));
      m += 1.0;
      ss[i*kk+j]=c1*sin(m*pi*(x-0.5));
    }
}


/*  Mass term, this diverges as beta goes to zero. */

static void mass22(element *cc, element *ss, int top, 
		   element beta, element offset){
  int i,j;
  element pre,pi,m,n;
  doublereal b2=2.0*beta;
  if(beta<10.0*DBL_EPSILON){
    for(i=0; i<top*top; i++)
	ss[i]=cc[i]=0.0;
  } else {
    pi=4.0*atan(1.0);
    /*  There is an extra factor of 0.5 to make definition standard */
    pre=0.125*sqrt(pi)*MYGAMMA(&b2);         /* 0.25 -> 0.125 */
    for(i=0; i<top; i++)
      for(j=0; j<=i; j++){
	m=2.0*(element) i+offset;
	n=2.0*(element) j+offset;
	cc[j*top+i]=cc[i*top+j]=pre*(bj(2*beta-0.5,pi*(m-n))+
		      bj(2*beta-0.5,pi*(m+n)));
	m += 1.0; n += 1.0;
	ss[j*top+i]=ss[i*top+j]=pre*(bj(2*beta-0.5,pi*(m-n))-
		      bj(2*beta-0.5,pi*(m+n)));
      }
  }
}

/*  lambda interaction */

static void lambda22(element *cc, element *ss, int top, 
		     element beta, element offset){
  int i,j;
  element pre,pi,m,n;
  doublereal b2=beta+0.5;
  pi=4.0*atan(1.0);
  /*  There is an extra factor of 1/(4 pi) to make definition standard */
  pre=0.125*pow(MYGAMMA(&b2),2);           /* 0.25 -> 0.125 */
  for(i=0; i<top; i++)
    for(j=0; j<=i; j++){
      m=2.0*(element) i+offset; n=2.0*(element) j+offset;
      cc[j*top+i]=cc[i*top+j]=pre*bj(beta,pi*m)*bj(beta,pi*n);
      ss[j*top+i]=ss[i*top+j]=0.0;
    }
}

/*  instantaneous annihilation interaction */

static void annihilate22(element *cc, element *ss, int top, 
			 element beta, element offset){
  int i,j;
  element pre,pi,m,n;
  doublereal b2=beta+0.5;
  pi=4.0*atan(1.0);
  /*  There is an extra factor of -1/(8 pi) to make definition standard */
  pre = -0.0625*pow(pi*MYGAMMA(&b2),2);  /* 0.125 -> 0.0625 */
  for(i=0; i<top; i++)
    for(j=0; j<=i; j++){
      cc[j*top+i]=cc[i*top+j]=0.0;
      m=2.0*(element) i+1.0+offset; n=2.0*(element) j+1.0+offset;
      ss[j*top+i]=ss[i*top+j]=
	pre*m*n*bj(beta+1.0,pi*m)*bj(beta+1.0,pi*n);
    }
}

/*  Matrix elements of the instantaneous interaction for 
    the sine and cosine bases.                           
    This is calculated as a sum of Gauss-Jacobi numerical integration
    and an excactly solved integral.    
    new version twice as fast as the older version shown later
    since there is less trig computed.  */

#define MAXGAUSS (4*TOP+8) /*Max number of gauss points for integration*/


static void instant22(element *cc, element *ss, int top, 
		      element beta, element offset){
  int i,j,rpc,rmc,row,col;
  element av,bv,au,bu,cau1,cau2,sav1,sav2,cbu,sbu,cbv,sbv,temp;
  element pi,u,v,a1,b1,vv;
  integer kind=5,nn,kpts=0;
  doublereal alpha=beta-0.5,beta2=beta+0.5,*lims=NULL,
    w[MAXGAUSS],t[MAXGAUSS],work[MAXGAUSS];
#ifdef USE_MPI /* here we create a local context for communication */
  double sendercc[(TOP+1)/2],senderss[(TOP+1)/2];
  int nslaves,myrank;
  MPI_Comm_size(IMPROVED_MPI_COMM,&nslaves);
  MPI_Comm_rank(IMPROVED_MPI_COMM,&myrank);
#if 0 /* debugging print */
  printf("Proc %i of %i starting instant 22\n",myrank,nslaves);
#endif
#endif
  pi=4.0*atan(1.0);

  /* rpc is the row number plus column number */
  for(rpc=0; rpc<2*top-1; rpc++)
#ifdef USE_MPI
    if(rpc%nslaves==myrank)
#endif
  {
    /* find weights and rescale interval */
    /* nn is the number of points and weights to use
       nn=10+2*rpc is equal to m+n+10 for cosine
       and m+n+8 for sine wavefunctions.  This should give about
       11 digits of accuracy.                */
    nn=10+2*rpc;
    if(nn>MAXGAUSS){
       fprintf(stderr,__FILE__":  Too many gauss points in instant22, exiting\n");
       fprintf(stderr,__FILE__":  top=%i rpc=%i MAXGAUSS=%i\n",top,rpc,MAXGAUSS);
      return;
    }

    /* calculate gaussian weights and points */
    GAUSSQ(&kind,&nn,&alpha,&alpha,&kpts,lims,work,t,w);
    for(i=0; i<nn; i++){
      t[i]=0.5*(1.0+t[i]);
      w[i]/=pow(2.0,2.0*alpha+1);
    }

    /* start with constant part */
    /* rmc is the row number minus column number */
    for(rmc=rpc%2; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2){
      row=(rpc+rmc)/2; col=(rpc-rmc)/2;
      /*  First, do cosine states  */
      a1=pi*((element) rpc+offset); 
      b1=pi*(element) rmc;
      cc[row*top+col] = -0.5*pi*pow(MYGAMMA(&beta2),2)*(
	 bj(beta,b1)*((1.0+4.0*beta)*bj(beta,b1)-2.0*b1*b1*bj(beta+1.0,b1))+
	 bj(beta,a1)*((1.0+4.0*beta)*bj(beta,a1)-2.0*a1*a1*bj(beta+1.0,a1)));
      /* second to sine states */
      a1 += pi;
      ss[row*top+col]= -0.5*pi*pow(MYGAMMA(&beta2),2)*(
	 bj(beta,b1)*((1.0+4.0*beta)*bj(beta,b1)-2.0*b1*b1*bj(beta+1.0,b1))-
	 bj(beta,a1)*((1.0+4.0*beta)*bj(beta,a1)-2.0*a1*a1*bj(beta+1.0,a1)));
#if 0 /* debugging print */
      if(rpc<2)printf("   mygamma %f bj %f bj %f\n",
		       MYGAMMA(&beta2),bj(beta+1.0,b1),bj(beta+1.0,a1));
      if(rpc<2)printf("   rpc=%i cc=%f ss=%f\n",
		       rpc,cc[row*top+col],ss[row*top+col]);
#endif
    }

    /* now do numerical integration */
    /* The loops are done in a funny order so that we don't spend
       all our time calculating trig functions */
    /*  First, do cosine states  */
    for(i=0; i<nn; i++){
      /* diagonal part of integration */
      u=2.0*t[i]-1.0;
      u=2.0*t[i]-1.0;
      vv = -0.25*w[i]*w[i]*(1.0-u*u);
      /* first do cosines */
      a1=pi*((element) rpc+offset); 
      cau1=vv*cos(a1*u);
      sav1=vv*a1*a1;
      /* then sines */
      a1+=pi; 
      cau2=vv*cos(a1*u);
      sav2=vv*a1*a1;
      au=2.0*pow(sin(pi*u),2); bu=sin(2.0*pi*u);
      rmc=rpc%2;
      cbu=cos(pi*rmc*u); sbu=sin(pi*rmc*u);
      for(; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2){
	row=(rpc+rmc)/2; col=(rpc-rmc)/2;
	b1=pi*(element) rmc;
	temp=b1*b1;
	cc[row*top+col] += cbu*sav1+cau1*temp;
	ss[row*top+col] += cbu*sav2-cau2*temp;
	temp=cbu;
	cbu=cbu-(au*cbu+bu*sbu);
	sbu=sbu-(au*sbu-bu*temp);
      } 
      /* non-diagonal part of integration (this should be what limits
	 speed).  Using symmetries, we only do 1/4 of the area 
	 of the square.  We also use recursion for the trig expressions. */
      for(j=0; j<i && i+j<nn; j++){
	u=t[i]+t[j]-1.0;  v=t[i]-t[j];
	/* include an extra factor of 1/2 for terms on other diagonal */
	if(i+j+1==nn)
	  vv = -2.0*w[i]*w[j]*(1.0-u*u)/(v*v);
	else
	  vv = -4.0*w[i]*w[j]*(1.0-u*u)/(v*v);
	/* first do cosines */
	a1=pi*((element) rpc+offset); 
	cau1=vv*cos(a1*u);
	sav1=vv*pow(sin(0.5*a1*v),2);
	/* then sines */
	a1+=pi;
	cau2=vv*cos(a1*u);
	sav2=vv*pow(sin(0.5*a1*v),2);
	au=2.0*pow(sin(pi*u),2); bu=sin(2.0*pi*u);
        av=2.0*pow(sin(0.5*pi*v),2); bv=sin(pi*v);
	rmc=rpc%2;
	cbu=cos(pi*rmc*u); sbu=sin(pi*rmc*u);
	cbv=cos(0.5*pi*rmc*v); sbv=sin(0.5*pi*rmc*v);
	for(; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2){
	  row=(rpc+rmc)/2; col=(rpc-rmc)/2;
	  temp=pow(sbv,2);
	  cc[row*top+col] += cbu*sav1+cau1*temp;
	  ss[row*top+col] += cbu*sav2-cau2*temp;
	  temp=cbu;
	  cbu=cbu-(au*cbu+bu*sbu);
	  sbu=sbu-(au*sbu-bu*temp);
	  temp=cbv;
	  cbv=cbv-(av*cbv+bv*sbv);
	  sbv=sbv-(av*sbv-bv*temp);
	}
      }
    }

    /* divide by -8.0 pi to give agreement with the notes */
    for(rmc=rpc%2; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2){
      row=(rpc+rmc)/2; col=(rpc-rmc)/2;
      cc[col*top+row]=cc[row*top+col] /= -16.0*pi;  /* 8->16 */
      ss[col*top+row]=ss[row*top+col] /= -16.0*pi;  /* 8->16 */
#if 0 /* debugging print */
      if(row<5 && col<5)
	printf(" row=%i col=%i cc=%.16f ss=%.16f\n",
	       row,col,cc[row*top+col],ss[row*top+col]);
#endif
    }
  }
#ifdef USE_MPI
#if 0  /* debugging print */
  printf("Proc %i of %i starting to collect from others\n",myrank,nslaves);
#endif 
  for(rpc=0; rpc<2*top-1; rpc++){
    /* we need to find i even if no data is broadcast */
    for(i=0, rmc=rpc%2; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2, i++)
      if(rpc%nslaves==myrank){
	row=(rpc+rmc)/2; col=(rpc-rmc)/2;
	sendercc[i]=cc[row*top+col];
	senderss[i]=ss[row*top+col];
      } 
    if(i>(TOP+1)/2){
       fprintf(stderr,__FILE__":  i=%i is too large in instant22\n",i);
      goto error;
    }  
    MPI_Bcast(sendercc,i,MPI_DOUBLE,rpc%nslaves,IMPROVED_MPI_COMM);
    MPI_Bcast(senderss,i,MPI_DOUBLE,rpc%nslaves,IMPROVED_MPI_COMM);
    if(rpc%nslaves!=myrank) 
      for(i=0, rmc=rpc%2; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2, i++){
	row=(rpc+rmc)/2; col=(rpc-rmc)/2;
	cc[row*top+col]=cc[col*top+row]=sendercc[i];
	ss[row*top+col]=ss[col*top+row]=senderss[i];
      }
  }
#endif
#if 0 /* print out some elements of ss and cc */
  for(i=0; i<5 && i<top; i++)
    for(j=0; j<5 && j<top; j++)
      printf("   %i %i %f %f\n",i,j,cc[i*top+j],ss[i*top+j]);
  printf("\n");
#endif
  return;  /* normal return */
  
#ifdef USE_MPI
 error: 
  fprintf(stderr,__FILE__":  Error in instant22, exiting\n");
  return;
#endif
}


#if 0  /* old version -slower (but gives identical results) */

static void instant22(element *cc, element *ss, int top, 
		      element beta, element offset){
  int i,j,rpc,rmc,row,col;
  element m,n;
  element pi,xc,u,v,a1,b1,vv[MAXGAUSS][MAXGAUSS/2];
  integer kind=5,nn,kpts=0;
  doublereal alpha=beta-0.5,beta2=beta+0.5,*lims=NULL,
    w[MAXGAUSS],t[MAXGAUSS],vvdiag[MAXGAUSS];
  pi=4.0*atan(1.0);

  /* rpc is the row number plus column number */
  for(rpc=0; rpc<2*top-1; rpc++){
    /* find weights and rescale interval */
    /* nn is the number of points and weights to use
       the following choice is equal to m+n+10 for cosine
       and m+n+8 for sine wavefunctions.  This should give about
       11 digits of accuracy.                */
    nn=10+2*rpc;
    if(nn>MAXGAUSS){
       fprintf(stderr,__FILE__":  Too many gauss points in instant22, exiting\n");
       fprintf(stderr,__FILE__":  top=%i rpc=%i MAXGAUSS=%i\n",top,rpc,MAXGAUSS);
      return;
    }
    /* calculate gaussian weights and points */
    /* vvdiag is used as a working array     */
    GAUSSQ(&kind,&nn,&alpha,&alpha,&kpts,lims,vvdiag,t,w);
    for(i=0; i<nn; i++){
      t[i]=0.5*(1.0+t[i]);
      w[i]/=pow(2.0,2.0*alpha+1);
    }
    /* calculate product of weights and interaction 
       (without trig factors) */
    for(i=0; i<nn; i++){
      u=2.0*t[i]-1.0;
      vvdiag[i] = -0.25*w[i]*w[i]*(1.0-u*u);
      for(j=0; j<i && i+j<nn; j++){
	u=t[i]+t[j]-1.0;
	v=t[i]-t[j];
	/* include an extra factor of 1/2 for terms on other diagonal */
	if(i+j+1==nn)
	  vv[i][j] = -2.0*w[i]*w[j]*(1.0-u*u)/(v*v);
	else
	  vv[i][j] = -4.0*w[i]*w[j]*(1.0-u*u)/(v*v);
      }
    }
    /* rmc is the row number minus column number */
    for(rmc=rpc%2; rmc<=rpc && rmc<=2*top-2-rpc; rmc+=2){
      row=(rpc+rmc)/2; col=(rpc-rmc)/2;
      /*  First, do cosine states  */
      m=2*(element) row+offset; n=2*(element) col+offset;
      a1=0.5*pi*(m+n); 
      b1=0.5*pi*(m-n);
      xc= -0.5*pi*pow(MYGAMMA(&beta2),2)*(
	 bj(beta,b1)*((1.0+4.0*beta)*bj(beta,b1)-2.0*b1*b1*bj(beta+1.0,b1))+
	 bj(beta,a1)*((1.0+4.0*beta)*bj(beta,a1)-2.0*a1*a1*bj(beta+1.0,a1)));
      /* numerical integration using gauss-Jacobi */
      for(i=0; i<nn; i++){
	/* diagonal part of integration */
	u=2.0*t[i]-1.0;
	xc += vvdiag[i]*(cos(b1*u)*a1*a1+cos(a1*u)*b1*b1); 
  /* non-diagonal part of integration (this should be what limits
    speed).  Using symmetries, we only do 1/4 of the area 
    of the square */
	for(j=0; j<i && i+j<nn; j++){
	  u=t[i]+t[j]-1.0;  v=t[i]-t[j];
	  xc += vv[i][j]*(cos(b1*u)*pow(sin(0.5*a1*v),2)+
			  cos(a1*u)*pow(sin(0.5*b1*v),2));
	}
      }
      /* divide by -8.0 pi to give agreement with the notes */
      cc[row*top+col]=cc[col*top+row]=-xc/(16.0*pi);  /* 8->16 */
      /* second to sine states */
      m += 1.0; n += 1.0;
      a1=0.5*pi*(m+n); b1=0.5*pi*(m-n);
      xc= -0.5*pi*pow(MYGAMMA(&beta2),2)*(
	 bj(beta,b1)*((1.0+4.0*beta)*bj(beta,b1)-2.0*b1*b1*bj(beta+1.0,b1))-
	 bj(beta,a1)*((1.0+4.0*beta)*bj(beta,a1)-2.0*a1*a1*bj(beta+1.0,a1)));
      /* numerical integration using gauss-Jacobi */
      for(i=0; i<nn; i++){
	/* diagonal part of integration */
	u=2.0*t[i]-1.0;
	xc += vvdiag[i]*(cos(b1*u)*a1*a1-cos(a1*u)*b1*b1); 
  /* non-diagonal part of integration (this should be what limits
    speed).  Using symmetries, we only do 1/4 of the area 
    of the square */
	for(j=0; j<i && i+j<nn; j++){
	  u=t[i]+t[j]-1.0; v=t[i]-t[j];
	  xc += vv[i][j]*(cos(b1*u)*pow(sin(0.5*a1*v),2)-
			  cos(a1*u)*pow(sin(0.5*b1*v),2));
	}
      }
      /* divide by -8.0 pi to give agreement with the notes */
      ss[row*top+col]=ss[col*top+row] = -xc/(16.0*pi);  /* 8->16 */
#if 1==0 /* debugging print */
      printf(" row=%i col=%i cc=%.16f ss=%.16f\n",
	     row,col,cc[row*top+col],ss[row*top+col]);
#endif
    }
  }
}

#endif

/* Calculate integrals for 1->3 particle interactions
   All momenta are the DLCQ values minus 0.5
   which=1 inner product 
   which=2 lambda interaction 
           (without the 1/(4 pi sqrt(K inner product))
   which=3 instantaneous interaction 
           (without the 1/(8 pi sqrt(K inner product))

   Unfortunately, the integration subroutine only
   uses relative errors.  Thus I have to do some gymnastics
   to get a consistant relative answer for the final
   answer.
   K=20 took 100 seconds on whiskey for tol=1.0e-6 (relative)
   K=20 took  61 seconds on whiskey for tol=1.0e-5 (relative)  
   
   This routine returns the error code of the integration routine.
   ierr=2 indicates error for DGAUS8.                      
*/

element inner13(int which, int kk,int k1,int k2,element b,integer *ierr){
  /* tol is the absolute error for the improved matrix elements */
  doublereal lower,upper,y,yt=0,tol=1.0e-10;
  f.kk=0.5+(element) kk;
  f.k1=0.5+(element) k1;
  f.k2=0.5+(element) k2;
  f.b=b;
  f.which=which;
  /* error flag, successful completion is ierr=0 */
#if 0 /* debug print */
  printf("starting inner13:  which=%i kk=%i k1=%i k2=%i\n",
	 which,kk,k1,k2);
#endif
#if 0  /*  This makes the above tol an absolute error */
  /*else tol is the relative error in the final improved
    value. */
  f.x=f.k1;
  tol/=f33(&f.k2);
#endif
  if(which==3 && k1==k2) return 0.0;
  /* Do K=3/2 in terms of K=5/2 */
  if(kk==1 && which==1)
      return 3.0*pow(3.0/5.0,2.0-6.0*b)*inner13(which,2,0,0,b,ierr);
  if(kk==1 && which==2)
      return 3.0*pow(3.0/5.0,0.5-3.0*b)*inner13(which,2,0,0,b,ierr);
  /* Do either of two lower corners */
  if((k1==0 && k2==0)||(k1+k2+1==kk && k2==0)){
    lower=0.0;
    upper=5.0/4.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F1_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr);
    yt+=y;
    lower=5.0/4.0;
    upper=5.0/3.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F1_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr);
    yt+=y;
    /* Do either lower or right edge */
  } else if((k2==0 || k1+k2+1==kk) && k1>0 && k1<kk-1){
#if 0 /* debugging print */
    printf("starting edge k1=%i/2, k2=%i/2, which=%i\n",
	   2*k1+1,2*k2+1,which);
#endif
    lower=f.k1-2.0/3.0;
    upper=f.k1-1.0/3.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F2_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr);
    yt+=y;
    lower=f.k1-1.0/3.0;
    upper=f.k1-1.0/4.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F2_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr);
    yt+=y;
    lower=f.k1-1.0/4.0;
    upper=f.k1+1.0/3.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F2_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr);
    yt+=y;
    lower=f.k1+1.0/3.0;
    upper=f.k1+3.0/4.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F2_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr);
    yt+=y;
    /* do interior region */
  } else if(k1>0 && k2>0 && k1+k2+1<kk){
    lower=f.k1-2.0/3.0; 
    upper=f.k1-1.0/3.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F3_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr); 
    yt+=y;
    lower=f.k1-1.0/3.0; 
    upper=f.k1+1.0/3.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F3_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr); 
    yt+=y;
    lower=f.k1+1.0/3.0; 
    upper=f.k1+2.0/3.0;
    f.tol=tol/(upper-lower);
    DGAUS8(F3_TO_FORTRAN,&lower,&upper,&f.tol,&y,ierr); 
    yt+=y;
  } else{
    fprintf(stderr,__FILE__":  incorrect input in inner13\n");
    fprintf(stderr,"           kk=%i k1=%i k2=%i\n",kk,k1,k2);
    goto error;
  }
  if(abs(*ierr)!=1)goto error;
  *ierr=0; /* no error */
  if(which<1 || which>3)goto error;

  return yt;  /* normal return */

error:
     fprintf(stderr,__FILE__":  Error in DGAUS8 in inner13, ierr=%li, which=%i\n",
	    *ierr,which);
    *ierr=2;
    return 0.0;  
}

#ifdef __GNUC__ /* assume that g77 is being used, doesn't have RECURSIVE */
#define DGAUS8_SUB DGAUS8_COPY
#else /* otherwise, call DGAUS8 recursively */
#define DGAUS8_SUB DGAUS8
#endif

/* Numerical integration function for corners 
   for which=1 and beta small, the integrand of
   the z integration disappears. */

doublereal f1(doublereal *u){
  integer ierr;
  doublereal lower=0.0,upper,y,tol;
  f.u=*u;
  y=2.0-5.0/(2.0* *u);
  if(y>lower)lower=y;
  if(f.which==3){
    upper=1;
    y=5.0/(2.0* *u)-1.0;
    if(y<upper)upper=y;
  } else
    upper=0.5;
  if(f.which==1 && f.b<10.0*DBL_EPSILON)
    y=0.0;
  else {
    tol=f.tol/(upper-lower);
    if(f.which==1) tol/=4.0*f.b;
    DGAUS8_SUB(F11_TO_FORTRAN,&lower,&upper,&tol,&y,&ierr);
    if(ierr==2) fprintf(stderr,__FILE__":  Error=%li in f1 DGAUSS8\n",ierr);
  }
  if(f.which==1)
    return 2.0* *u*(y+(pow(upper,2.0*f.b+1.0)-pow(lower,2.0*f.b+1.0))/
		    ((1.0+2.0*f.b)*pow(f.kk*f.kk*(f.kk-*u),2.0*f.b)));
  else if(f.which==2)
    return 2.0*(y+(pow(upper,f.b+0.5)-pow(lower,f.b+0.5))/
		((f.b+0.5)*pow(f.kk,2.0*f.b)*pow(f.kk-*u,f.b+0.5)));
  else
    return y+
      (f.kk-2.0* *u)*(pow(upper,f.b+0.5)-pow(lower,f.b+0.5))/
       ((f.b+0.5)*pow(f.kk,2.0*f.b+1.0)*pow(f.kk-*u,f.b+0.5))+
       (f.kk+ *u)*(pow(1.0-lower,f.b+0.5)-pow(1.0-upper,f.b+0.5))/
       ((f.b+0.5)*pow(f.kk,2.0*f.b)*pow(f.kk-*u,f.b+1.5));  
}

doublereal f11(doublereal *z){
  if(f.which==1)
    return pow((f.kk-f.u)*(1.0-*z)* *z,2.0*f.b)/
      pow((f.kk-f.u* *z)*(f.kk-f.u+f.u * *z),4.0*f.b)-
      pow(*z/(f.kk*f.kk*(f.kk-f.u)),2.0*f.b);
  else if(f.which==2)
    return pow((f.kk-f.u)*(1.0-*z)* *z,f.b-0.5)/
      pow((f.kk-f.u* *z)*(f.kk-f.u+f.u * *z),2.0*f.b)-
      pow(*z,f.b-0.5)/(pow(f.kk,2.0*f.b)*pow(f.kk-f.u,f.b+0.5));
  else
    return (f.kk+f.u* *z)*(f.kk-2.0*f.u+f.u* *z)*
      pow((f.kk-f.u)*(1.0-*z)* *z,f.b-0.5)/
      (pow(f.kk-f.u* *z,2.0*f.b+2.0)*pow(f.kk-f.u+f.u* *z,2.0*f.b))-
      (f.kk-2.0*f.u)*pow(*z,f.b-0.5)/
      (pow(f.kk,2.0*f.b+1.0)*pow(f.kk-f.u,f.b+0.5))-
      (f.kk+f.u)*pow(1.0-*z,f.b-0.5)/
      (pow(f.kk,2.0*f.b)*pow(f.kk-f.u,f.b+1.5));
}

/* integration routines for the lower and right edges

   Since we add and subtract the function at the 
   edge, the error estimation is more complicated.
   The size of the numerically evaluated integral
   to the full integral is estimated using the 
   next order perturbation in y.                  */

doublereal f2(doublereal *x){
  integer ierr;
  doublereal lower=0,upper,yt,tol;
  f.x=*x;
  yt=2.0*f.k1-0.5-2.0* *x;
  if(yt>lower)lower=yt;
  upper=1.5-f.k1+ *x;
  yt=0.5*f.k1+1.0-0.5* *x;
  if(upper>yt)upper=yt;
  yt=2.0*f.k1+1.5-2.0* *x;
  if(upper>yt)upper=yt;
  if(f.which==1 && f.b<10.0*DBL_EPSILON)
    yt=0.0;
  else {
    if(f.x<f.kk-f.x)
      tol=f.x;
    else
      tol=f.kk-f.x;
    tol*=f.tol/(upper-lower);
    if(f.which==1) tol/=4.0*f.b;
    DGAUS8_SUB(F22_TO_FORTRAN,&lower,&upper,&tol,&yt,&ierr);
    if(ierr==2) fprintf(stderr,__FILE__":  Error=%li in f2 DGAUSS8\n",ierr);
  }
  if(f.which==1)
    return yt+(pow(upper,2.0*f.b+1.0)-pow(lower,2.0*f.b+1.0))/
      ((1.0+2.0*f.b)*pow(f.kk*f.kk* *x*(f.kk-*x),2.0*f.b));
  else if(f.which==2)
    return yt+(pow(upper,f.b+0.5)-pow(lower,f.b+0.5))/
      ((f.b+0.5)*pow(f.kk,2.0*f.b)*pow(*x*(f.kk-*x),f.b+0.5));
  else if(fabs(f.k2-0.5)<0.1)
      return yt+(2.0*f.kk-*x)*(pow(upper,f.b+0.5)-pow(lower,f.b+0.5))/
	((f.b+0.5)*pow(f.kk,2.0*f.b)*pow(*x,f.b+1.5)*pow(f.kk-*x,f.b+0.5));
    else
      return yt+(2.0* *x-f.kk)*(pow(upper,f.b+0.5)-pow(lower,f.b+0.5))/
	((f.b+0.5)*pow(f.kk,2.0*f.b+1.0)*pow(*x*(f.kk-*x),f.b+0.5));
}

/*  In practice, the following is the largest 
    time consumer for K=20    */

doublereal f22(doublereal *y){
  doublereal a,b,c;
  if(f.which==1){
    a=(f.x+*y)*(f.kk-f.x)*(f.kk-*y);
    return pow((f.kk-f.x-*y)*f.x* *y/(a*a),2.0*f.b)-
      pow(*y/(f.kk*f.kk*f.x*(f.kk-f.x)),2.0*f.b);
  } else if(f.which==2){
    a=(f.x+*y)*(f.kk-f.x)*(f.kk-*y);
    b=f.kk*f.x*(f.kk-f.x);
    return pow((f.kk-f.x-*y)*f.x* *y/(a*a),f.b-0.5)/a-
      pow(*y/(f.kk*b),f.b-0.5)/b;
  } else if(fabs(f.k2-0.5)<0.1){
    a=f.x+*y;
    b=a*(f.kk-f.x)*(f.kk-*y);
    c=f.kk*f.x*(f.kk-f.x);
    return (f.x-*y)*(2.0*f.kk-f.x-*y)*
      pow((f.kk-f.x-*y)*f.x* *y/(b*b),f.b-0.5)/(a*a*b)-
      (2.0*f.kk-f.x)*pow(*y/(f.kk*c),f.b-0.5)/(f.x*c);
  } else {
    a=f.kk-*y;
    b=a*(f.x+*y)*(f.kk-f.x);
    c=f.kk*f.kk*f.x*(f.kk-f.x);
    return (2.0*f.x-f.kk+*y)*(f.kk+*y)*
      pow((f.kk-f.x-*y)*f.x* *y/(b*b),f.b-0.5)/(a*a*b)-
      (2.0*f.x-f.kk)*pow(*y/c,f.b-0.5)/c;
  }
}

doublereal f3(doublereal *x){
  integer ierr;
  doublereal lower,upper,yt,tol;
  f.x=*x;
  lower=2.0*f.k1+f.k2-1.0-2.0* *x;
  yt=0.5*f.k1+f.k2-0.5-0.5* *x;
  if(yt>lower)lower=yt;
  yt=*x-1.0-f.k1+f.k2;
  if(yt>lower)lower=yt;
  upper=2.0*f.k1+f.k2+1.0-2.0* *x;
  yt=0.5*f.k1+f.k2+0.5-0.5* *x;
  if(upper>yt)upper=yt;
  yt=*x+1.0-f.k1+f.k2;
  if(upper>yt)upper=yt;
  tol=f.tol/(upper-lower);
  DGAUS8_SUB(F33_TO_FORTRAN,&lower,&upper,&tol,&yt,&ierr);
  if(ierr==2) fprintf(stderr,__FILE__":  Error=%li in f3 DGAUSS8\n",ierr);
  return yt;
}

/* This should be the dominant time consumer
   for large K for the 1->3 particle interactions */

doublereal f33(doublereal *y){
  doublereal a,b;
  if(f.which==1){
    a=(f.x+*y)*(f.kk-f.x)*(f.kk-*y);
    return pow((f.kk-f.x-*y)*f.x* *y/(a*a),2.0*f.b);
  } else if(f.which==2){
    a=(f.x+*y)*(f.kk-f.x)*(f.kk-*y);
    return pow(*y*f.x*(f.kk-f.x-*y)/(a*a),f.b-0.5)/a;
  } else {
    a=f.x+*y;
    b=a*(f.kk-f.x)*(f.kk-*y);
    return (f.x-*y)*(2.0*f.kk-f.x-*y)*
      pow(*y*f.x*(f.kk-f.x-*y)/(b*b),f.b-0.5)/
      (a*a*b);
  }
}

#endif  /* functions for IMPROVED=3 */

/* This prints out the memory used by fun22 and fun13 */

void print_interaction_memory(FILE *fp){
  fprintf(fp,"MM 2->2 and 1->3 interactions use " SIZE_T_PRINT 
	  " and " SIZE_T_PRINT " bytes\n",
	      fun22_memory,fun13_memory);
  fprintf(fp,"MM   with pointer static memory use %i bytes\n",
	      sizeof(fun22)+sizeof(fun13));
}

/* 
   free memory allocated for fun22 and fun13. 
*/

void free_improved(void){
  int i,j,kkdown;
  
  /* kkdown is the sum of the momenta, *each* rounded down */
  for(kkdown=0; kkdown<=global_kt/2-(GLUE_PERIODIC?0:1); kkdown++)
    /* j is one incoming momentum (rounded down) */
    for(j=0; j<=kkdown; j++){
      free(fun22[kkdown][j]);
      fun22[kkdown][j]=NULL;
    }
  fun22_memory=0;
  last_kt=0;

  /* i is the momentum (rounded down) of the single incoming particle */
  for(i=(GLUE_PERIODIC?0:1); i<=(global_kt-(GLUE_PERIODIC?0:1))/2; i++)
    /* j is the momentum (rounded down) of an outgoing particle */
    for(j=0; j<=(GLUE_PERIODIC?i:i-1); j++){
      free(fun13[i][j]);
      fun13[i][j]=NULL;
    }
  fun13_memory=0;
  last_beta=-1.0;
  last_kt=0;
}