//34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
// Subroutines for numerical integration of one loop Feynman diagrams.
// D.E. Soper
// 25 September 2006
//
#include "Nphoton.h"
#include <cmath>
#include <iostream>
#include <cassert>
using namespace std;

//==================================================================================================
// Implementation code for class DESxpoint.
// Combination of point.neu and point.rho tested by integrating f = const. with
// psoft = 0.9, cparam = 6, 5, 4, 3, 2, 1, 0.5 to accuracy 0.002,  20 November 2005, DES.
// 9 March 2006, DES.
//
// Constructor. Note that S needs to be specified in DESxpoint::newS(DES_Smatrix sIN).
DESxpoint::DESxpoint(int seedIN, int nvertsIN, double xcutoffIN, double alphaIN[])
 : r(seedIN)   // Construct random number generator instance r.
{
seed   = seedIN ;
nverts = nvertsIN ;
nvertsp1 = nverts + 1;
xcutoff = xcutoffIN;
int i,n;
double alphasum;
//
alphasum = 0.0;
for(i = 1; i <= 11; i++) 
  {
  alpha[i] = alphaIN[i];
  alphasum = alphasum + alphaIN[i];
  }
assert(abs(alphasum - 1.0) <= 1.0e-10);
// The alpha[i] are
// 1: Probability for uniform distribution.
// 2: Probability to choose points with one of "soft-n" distributions.
// 3: Probability to choose points with one of the modified collinear distributions.
// 4: Probability to choose points with one of the modified soft distributions.
// 5: Probability to choose points with one of the "x0" collinear distritutions.
// 6: Probability to choose points with one of the "x^0 x^n x^np1" collinear distributions.
// 7: Probability to choose points with one of the collinear-soft distributions.
// 8: Probability to choose points with the double parton scattering distribution.
// 9: Probability to choose points with one of the colliner distribtions.
//10: Probability to choose points with small x^0, uniform x^i distribution.
//11: Probability to choose points with large x^0, uniform x^i distribution.
//
cparam = 0.5; // Controls power: eg. rho_n ~ 1/g_n(x)^{N-1-cparam} in soft-n distribution.
c10 = 0.5;    // Controls power for method 10; smaller c10 -> points more concentrated at small x^0.
c11 = 0.8;    // Controls power for method 11; smaller c11 -> points more concentrated at large x^0.
xsmax = 0.1;                  // Maximum xs in "x^0 x^n x^np1" collinear distributions.
for(n = 0; n <= nverts; n++)  // Parameters A_n for restriction g_n(x) < A_n in chosing points.
  {
  capA[n] = 1.0;
  }
}

// There is no default constructor.
//--------------------------------------------------------------------------------------------------
// Change matrix S.
void DESxpoint::newS(DES_Smatrix sIN)
{
int i,j,n,np1,naa,nbb,aa,bb,aap1,bbp1;
double det;
double min_s,temp,sproduct;
double wsoft[nverts];
bool signsOK;
double denom1,denom2,denom3,denom4,df;
double epsilondps = 0.3;
double epsiloncoll = 0.3;
//
s0 = sIN.getrtssq();        // Parameter with dimensions momentum^2, cancels from physical results.
m0sq = sIN.getm0sq();       // Squared mass parameter used to control UV convergence.
for(i = 0; i <= nverts; i++)
  {
  for(j = 0; j <= nverts; j++) s[i][j] = sIN(i,j);
  }
Ncollx0 = 0.0;
for(n = 1; n <= nverts; n++)  // Number of cases for the "x0" collinear distritutions.
  {
  if(s[cyclic(n-1,nverts)][cyclic(n+1,nverts)] > 0) Ncollx0 = Ncollx0 + 2.0;
  }
// Construct max values for variables in the collinear-soft schemes.
xAminusmax[0]  = 0.0;
xBminusmax[0]  = 0.0;
xnminusmax[0]  = 0.0;
xAplusmax[0]   = 0.0;
xBplusmax[0]   = 0.0;
xnp1plusmax[0] = 0.0;
for(n = 1; n <= nverts; n++)
  {
  xAminusmax[n]  = 0.5;
  xAplusmax[n]   = 0.5;
  for(i = 1; i <= nverts; i++)
    {
	if((i != cyclic(n-1,nverts)) && (i != n) && (i != cyclic(n+1,nverts)))
	  {
	  temp = 0.5*abs(s[i][n]/s[cyclic(n-1,nverts)][cyclic(n+1,nverts)]);
	  if(temp < xAminusmax[n]) xAminusmax[n] = temp;
	  }
	if((i != cyclic(n+2,nverts)) && (i != n) && (i != cyclic(n+1,nverts)))
	  {
	  temp = 0.5*abs(s[i][cyclic(n+1,nverts)]/s[n][cyclic(n+2,nverts)]);
	  if(temp < xAplusmax[n]) xAplusmax[n] = temp;
	  }
	}
  xBminusmax[n]  = 3.0;
  xBplusmax[n]   = 3.0;
  xnminusmax[n]  = 0.5;
  xnp1plusmax[n] = 0.5;
  }
// Construct coefficients a_ni and b_n.
b[0] = 0.0;                                 // These should never be called.
for(i = 0; i <= nverts; i++) a[0][i] = 0.0;  // These should never be called.
// Now for each n, find min of the |S_ni| and |S_{n-1}{n+1}| and use it to define a_ni and b_n.
wsoftsum = 0.0; // Initialization for sum of weights w_n^(s) = wsoft[n].
for(n = 1; n <= nverts; n++)
  {
  min_s = abs( s[cyclic(n-1,nverts)][cyclic(n+1,nverts)] );
  sproduct = min_s;
  for(i = n + 2; i <= n - 2 + nverts; i++)
    {
    temp = abs( s[n][cyclic(i,nverts)] );
    sproduct = sproduct*temp;
    if(temp < min_s) min_s = temp;
    }
  temp = m0sq;
  sproduct = sproduct*temp;
  if(temp < min_s) min_s = temp;
  lambda[n] = s0/min_s;          // Normalizing factor \lambda_n relating a_ni and b_n to |S_ij|
  // Calculate weights w_n^(s) for chosing index n for chosing a point x.
  wsoft[n] = pow(s0,nverts-1)*pow(capA[n]/lambda[n],cparam)/sproduct; 
  wsoftsum = wsoftsum + wsoft[n];
  b[n] = abs( s[cyclic(n-1,nverts)][cyclic(n+1,nverts)] )/min_s;
  a[n][0] = m0sq/min_s;          // This is used.
  for(i = 1; i <= nverts; i++) a[n][i] = abs(s[n][i])/min_s;
  a[n][cyclic(n-1,nverts)] = 0.0; // These should never be called.
  a[n][n] = 0.0;                  // These should never be called.
  a[n][cyclic(n+1,nverts)] = 0.0; // These should never be called.
  }
// Calculate partialprob[n] = sum_{n' = 1}^n  wsoft[n'] / sum_{n' = 1}^nverts  wsoft[n'].
partialprob[0] = 0.0;
for(n = 1; n <= nverts - 1; n++) partialprob[n] = partialprob[n-1] + wsoft[n]/wsoftsum;
partialprob[nverts] = 1.0;
// Generate information for double parton scattering method.
dps.n = 0;
if(nverts >= 6)
  {
  for(aa = 1; aa <= nverts - 3; aa++)
    {
    aap1 = aa + 1;
    for(bb = aa + 3; bb <= nverts && bb <= aa + nverts - 3; bb++)
      {
      bbp1 = bb + 1;
      if (bbp1 > nverts) bbp1 = bbp1 - nverts;
      signsOK = false;
      if(s[aap1][bb] > 0)
        {
        if(s[aa][bbp1] > 0 && s[aap1][bbp1] < 0 && s[aa][bb] < 0) signsOK = true;
        }
      else if (s[aap1][bb] < 0)
        {
        if(s[aa][bbp1] < 0 && s[aap1][bbp1] > 0 && s[aa][bb] > 0) signsOK = true;
        }
      if(signsOK)
        {
        assert(s[aap1][bb] > 0);
        det = s[aap1][bb]*s[aa][bbp1] - s[aap1][bbp1]*s[aa][bb];
        denom1 = s[aap1][bb] - s[aa][bb];
        denom2 = s[aa][bbp1] - s[aap1][bbp1];
        denom3 = s[aa][bbp1] - s[aa][bb];
        denom4 = s[aap1][bb] - s[aap1][bbp1];
        df = abs(det)/pow(denom1 + denom2,2);
        if(df < epsilondps)  // We have a good choice for aa and bb!
          {
          dps.n = dps.n + 1;
          dps.A = aa;
          dps.Ap1 = aap1;
          dps.B = bb;
          dps.Bp1 = bbp1;
          dps.df = df;
          dps.fA = denom4/(denom1 + denom2);
          dps.fB = denom2/(denom3 + denom4);
          dps.oneminusfA = 1.0 - dps.fA;
          dps.oneminusfB = 1.0 - dps.fB;
          assert(dps.n < 2);
          } // Close if(df < epsilondps).
        } // Close if (signsOK).
      }   // Close for(bb = aa + 4; bb <= aa - 4 + nverts; b++)
    }  // Close for (aa = 1; aa <= nverts - 1; aa++).
  } // Close if (nverts >= 6).
// Generate information for collinear singularity method.
collN = 0;
if(nverts >= 6)
  {
  for(n = 1; n <= nverts - 1; n++)
    {
    np1 = n + 1;
    for(naa = n + 3; naa <= n + nverts - 3; naa++)
      {
      for(nbb = naa + 1; nbb <= n + nverts - 2; nbb++)
        {
        aa = naa;
        if(aa > nverts) aa = aa - nverts;
        bb = nbb;
        if(bb > nverts) bb = bb - nverts;
        signsOK = false;
        if(s[n][aa] > 0)
          {
          if(s[np1][bb] > 0 && s[np1][aa] < 0 && s[n][bb] < 0) signsOK = true;
          }
        else if (s[n][aa] < 0)
          {
          if(s[np1][bb] < 0 && s[np1][aa] > 0 && s[n][bb] > 0) signsOK = true;
          }
        if(signsOK)
          {
          det = s[np1][aa]*s[n][bb] - s[np1][bb]*s[n][aa];
          denom1 = s[np1][aa] - s[n][aa];
          denom2 = s[n][bb] - s[np1][bb];
          df = abs(det)/denom1/denom2;
          if(df < epsiloncoll)  // We have a good choice for aa and bb!
            {
            collN = collN + 1;
            coll[collN].N = n;
            coll[collN].Np1 = np1;
            coll[collN].A = aa;
            coll[collN].B = bb;
            coll[collN].df = df;
            coll[collN].f = (s[np1][aa] - s[np1][bb])/(denom1 + denom2);
            } // Close if(df < epsiloncoll).
          }
        } // Close for(nbb = naa + 1; bb <= n + nverts - 2; nbb++)
      } // Close for(naa = n + 3; naa <= n + nverts - 3; naa++).
    } // Close for(n = 1; n <= nverts - 1; n++).
  }// Close if(nverts >= 6)
}
//--------------------------------------------------------------------------------------------------
// Generate new point x.
// Here x is an array x[0],x[1],...,x[nverts] with sum x[i] = 1.
void DESxpoint::neu(DESxarray& xOUT, bool& xOK, int& method)
{
static int i,j,n,nn,n1,n2,n3,nmax;
double rval[nverts];   // {r[0], r[1],..., r[nverts - 1]}
double yy[nverts + 1]; // {y[0], y[1],..., y[nverts]}
static double temp, rmethod;
static double omega,omegamax,omegamin,expomega;
static double xbar,sqrtxbar,ybar,xs,ys,yn1,rtxs;
static double sum;
static double r1,r2,rmax,rmin;
static bool minuschoice;
static double svalue;
static double xA,xB,y0;
static int nm1,np1,np2;
static double t,u;
static double oneminusxbar;
static double rbar,axbar,xbarpower;
static double tildeomega;
static double signA,signB,omegafactorA,omegafactorB;
static double sqrtxt;
static double xu,xt,oneminusxs;
static double xumax,deltau;
static double yA,yAp1,yB,yBp1;
static double tu,xa,rcoll,ycoll,dfsq,oneminusxt,xbarsq,acoll;
static double yN,yNp1;
static double w,onemx0;
static double lambdaR;
static double alphasum;
static const double twothirds = 0.6666666666666667;
static const double tiny = 1.0e-10;
//
// The alpha[i] are
// 1: Probability for uniform distribution.
// 2: Probability to choose points with one of "soft-n" distributions.
// 3: Probability to choose points with one of the modified collinear distributions.
// 4: Probability to choose points with one of the modified soft distributions.
// 5: Probability to choose points with one of the "x0" collinear distritutions.
// 6: Probability to choose points with one of the "x^0 x^n x^np1" collinear distributions.
// 7: Probability to choose points with one of the collinear-soft distributions.
// 8: Probability to choose points with the double parton scattering distribution.
// 9: Probability to choose points with one of the collinear distribtions.
//10: Probability to choose points with small x^0, uniform x^i distribution.
//11: Probability to choose points with large x^0, uniform x^i distribution.
//
for(n = 0; n <= nverts; n++) x[n] = 777.0; // Dummy value;
OK = false;
// Choose method.
rmethod = r.random();
method = 777; // dummy value
alphasum = 0.0;
for(i = 1; i <= 11; i++)
  {
  alphasum = alphasum + alpha[i];
  if(rmethod <= alphasum)
    {
    method = i;
    break;
    }
  }
if (method == 8 && dps.n == 0) method = 2; // Realocate to soft-n method if dpf method not needed.
if (method == 9 && collN == 0) method = 2; // Realocate to soft-n method if coll method not needed.
//
if(method == 1)
  {
  // First method, with uniform distribution
  for(i = 0; i < nverts; i++) rval[i] = r.random(); // Get nverts random numbers in (0,1).
  DESsort(rval,nverts);                              // Sort the random numbers.
  x[1] = rval[0];                // Construct x[1],..., x[nverts] as differences.
  for(i = 2; i <= nverts; i++)  
    {
    x[i] = rval[i-1] - rval[i-2];
    }
  x[0] = 1.0 - rval[nverts - 1];  // x[0] = 1 - sum x[i]
  OK = true;
  }
else if(method == 10)
  {
  // Well, this is the 10th method, but it is similar to the uniform distriubtion.
  for(i = 0; i < nverts - 1; i++) rval[i] = r.random(); // Get nverts - 1 random numbers in (0,1).
  DESsort(rval,nverts - 1);                             // Sort the random numbers.
  yy[2] = rval[0];                // Construct y[1],..., y[nverts] as differences.
  for(i = 3; i <= nverts; i++)  
    {
    yy[i] = rval[i-2] - rval[i-3];
    }
  yy[1] = 1.0 - rval[nverts - 2];
  w = 1.0 - pow(r.random(),1.0/c10);
  onemx0 = pow(w,1.0/nverts);
  x[0] = 1.0 - onemx0;
  for(i = 1; i <= nverts; i++)
    {
    x[i] = yy[i]*onemx0;
    }
  OK = true;
  // Density is c10*nverts!/(1-w)^{1-c10}.
  }
else if(method == 11)
  {
  // Well, this is the 11th method, but it is similar to the uniform distriubtion.
  for(i = 0; i < nverts - 1; i++) rval[i] = r.random(); // Get nverts - 1 random numbers in (0,1).
  DESsort(rval,nverts - 1);                             // Sort the random numbers.
  yy[2] = rval[0];                // Construct y[1],..., y[nverts] as differences.
  for(i = 3; i <= nverts; i++)  
    {
    yy[i] = rval[i-2] - rval[i-3];
    }
  yy[1] = 1.0 - rval[nverts - 2];
  onemx0 = pow(r.random(),1.0/c11);
  x[0] = 1.0 - onemx0;
  for(i = 1; i <= nverts; i++)
    {
    x[i] = yy[i]*onemx0;
    }
  OK = true;
  // Density is c11*(nverts-1)!/(onemx0)^{nverts-c11}.
  }
else if(method == 2)
  {
  // Second method, one of the soft-n distributions.
  // First choose n.
  temp = r.random(); 
  for(n = 1; partialprob[n] < temp; n++){}
  // This exits with n = first n with partialprob[n] >= temp.
  // Now chose the point;
  for(i = 0; i < nverts - 2; i++) rval[i] = r.random();  // Get nverts - 2 random numbers in (0,1).
  DESsort(rval,nverts - 2);                               // Sort the random numbers.
  yy[0] = rval[0];
  for(i = 1; i <= nverts - 3; i++)
    {
    nn = n + 1 + i;
    if(nn > nverts) nn = nn - nverts;  // Cyclic notation in {1, 2,..., nverts}.
    yy[nn] = rval[i] - rval[i-1];
    }
  ybar = 1.0 - rval[nverts - 3];
  //
  // Find x_s and use it to define xbar and the x_i, then omega
  xs = capA[n]*pow(r.random(),1.0/cparam);
  xbar = xs/b[n]*ybar;
  x[0] = xs/a[n][0]*yy[0];
  for(j = 1; j <= nverts - 3; j++)
    {
    i =  n + 1 + j;
    if(i > nverts) i = i - nverts;  // This covers i in {1,..., nverts} omitting n-1,n,n+1.
    x[i] = xs/a[n][i]*yy[i];
    }
  sqrtxbar = sqrt(xbar);
  omegamax = log(1.0/sqrtxbar);
  omega = (2.0*r.random() - 1.0)*omegamax;
  expomega = exp(omega);
  nn = n - 1;
  if(nn < 1) nn = nn + nverts;
  x[nn] = sqrtxbar*expomega;
  nn = n + 1;
  if(nn > nverts) nn = nn - nverts;
  x[nn] = sqrtxbar/expomega;
  // Finally, we need x[n].
  sum = 0.0;
  x[n] = 0.0; // So i = n is not included in sum below.
  for(i = 0; i <= nverts; i++) sum = sum + x[i];
  x[n] = 1.0 - sum;
  if(x[n] > 0.0) OK = true;
  else           OK = false;
  }
else if(method == 3)
  {
  // Third method, one of the modified collinear distributions.
  // x[0] + x[n] + x[n+1] is typically small.
  // First choose n.
  temp = r.random(); 
  for(n = 1; double(n)/double(nverts) < temp; n++){}
  // This exits with n = first n with n/nverts >= temp.
  for(i = 0; i < nverts - 2; i++) rval[i] = r.random();  // Get nverts - 2 random numbers in (0,1).
  DESsort(rval,nverts - 2);                               // Sort the random numbers.
  // Get x[n+2], x[n+3],..., x[n-1]  (in a cyclic notation).
  nn = n + 2;
  if(nn > nverts) nn = nn - nverts;
  x[nn] = rval[0];                     // Cyclic notation in {1, 2,..., nverts}.
  for(i = 1; i <= nverts - 3; i++)
    {
    nn = n + 2 + i;
    if(nn > nverts) nn = nn - nverts;  // Cyclic notation in {1, 2,..., nverts}.
    x[nn] = rval[i] - rval[i-1];
    }
  ybar = 1.0 - rval[nverts - 3]; // = 1 - (x[n+2] +  x[n+3] + ... + x[n-1])
  // Now we need x[0], x[n], x[n+1], which sum to ybar.
  r1 = r.random();
  r2 = r.random();
  if(r1 > r2)
    {
    rmax = r1;
    rmin = r2;
    }
  else
    {
    rmax = r2;
    rmin = r1;
    }
  x[0] = ybar*rmin;
  x[n] = ybar*(rmax - rmin);
  nn = n+1;
  if(nn > nverts) nn = nn - nverts;
  x[nn] = ybar*(1.0 - rmax);
  OK = true;
  }
else if(method == 4)
  {
  // Fourth method, one of the modified soft distributions.
  // x[n-1] + x[n] + x[n+1] is typically small.
  // First choose n.
  temp = r.random(); 
  for(n = 1; double(n)/double(nverts) < temp; n++){}
  // This exits with n = first n with n/nverts >= temp.
  for(i = 0; i < nverts - 2; i++) rval[i] = r.random(); // Get nverts - 2 random numbers in (0,1).
  DESsort(rval,nverts - 2);                              // Sort the random numbers.
  // Get x[0] and then x[n+2], x[n+3],..., x[n-2]  (in a cyclic notation).
  x[0] = rval[0];
  for(i = 1; i <= nverts - 3; i++)
    {
    nn = n + 1 + i;
    if(nn > nverts) nn = nn - nverts;  // Cyclic notation in {1, 2,..., nverts}.
    x[nn] = rval[i] - rval[i-1];
    }
  ybar = 1.0 - rval[nverts - 3]; // = 1 - x[0] - (x[n+2] +  x[n+3] + ... + x[n-2])
  // Now we need x[n-1], x[n], x[n+1], which sum to ybar.
  r1 = r.random();
  r2 = r.random();
  if(r1 > r2)
    {
    rmax = r1;
    rmin = r2;
    }
  else
    {
    rmax = r2;
    rmin = r1;
    }
  x[n] = ybar*rmin;
  nn = n - 1;
  if(nn < 1) nn = nn + nverts;
  x[nn] = ybar*(rmax - rmin);
  nn = n+1;
  if(nn > nverts) nn = nn - nverts;
  x[nn] = ybar*(1.0 - rmax);
  OK = true;
  }
  else if(method == 5)  
  {
  // Fifth method, one of the "x0" collinear distributions.
  // First we need to decide to choose points according to cs-minus[n] or cs-plus[n].
  temp = r.random(); 
  if(temp > 0.5) // We look for cs-minus choice. Now we have to choose n.
    {
    minuschoice = true;
	temp = r.random();
	i = 0;
	n = 1;
	while(true)
	  {
	  assert(n <= nverts);
	  if(s[cyclic(n-1,nverts)][cyclic(n+1,nverts)] > 0) i++; 
      // Allowed n choices have s[n-1][n+1] > 0.
	  if(double(2*i)/Ncollx0 > temp) // This is the n we want.
	    {
        n1 = n;
        n2 = cyclic(n-1,nverts);
        n3 = cyclic(n+1,nverts);
        nmax = n1;   // The greater of n1 and n2.
        break;       // We exit the while statement.
		}
	  n++;
	  }
	}
  else // We look for cs-plus choice. Now we have to choose n.
	{
	minuschoice = false;
	temp = r.random();
	i = 0;
	n = 1;
	while(true)
	  {
	  assert(n <= nverts);
	  if(s[cyclic(n,nverts)][cyclic(n+2,nverts)] > 0) i++;
      // Allowed n choices have s[n-1][n+1] > 0.
	  if(double(2*i)/Ncollx0 > temp) // This is the n we want.
	    {
        n1 = cyclic(n+1,nverts);
        n2 = cyclic(n+2,nverts);
        n3 = n;
        nmax = n2;  // The greater of n1 and n2.
        break;  // We exit the while statement.
        }
	  n++;
	  }
	} // close if(temp > 0.5) ... else ... .
  // Find x_s.
  xs = pow(r.random(),1.0/cparam);
  rtxs = sqrt(xs);
  // Find the y^j.
  for(i = 0; i < nverts - 3; i++) rval[i] = r.random(); // Get nverts - 3 random numbers in (0,1).
  DESsort(rval,nverts - 3);                              // Sort the random numbers.
  x[cyclic(nmax+1,nverts)] = xs*rval[0];      // Construct x[j] for j notin S as xs*differences.
  for(i = 2; i <= nverts - 3; i++)  
    {
    x[cyclic(nmax+i,nverts)] = xs*(rval[i-1] - rval[i-2]);
    }
  x[cyclic(nmax-2,nverts)] = xs*(1.0 - rval[nverts - 4]);  // x[n1 - 1] = xs - sum x[i]
  // Find x^0.
  x[0] = rtxs*((1.0 - xs + rtxs)/((1.0 - xs)*r.random() + rtxs) - 1.0);
  // Find x[n1] and x[n2].
  svalue = s[n2][n3];
  ys = svalue/(svalue  + 2.0*m0sq);
  omegamin = atan(- ys/rtxs);
  omegamax = atan((1.0 - ys)/rtxs);
  omega = (omegamax - omegamin)*r.random() + omegamin;
  yn1 = ys + rtxs*tan(omega);
  x[n1] = (1.0 - x[0] - xs)*yn1;
  x[n2] = (1.0 - x[0] - xs)*(1.0 - yn1);
  OK = true;
  }
else if(method == 6)
  {
  // Sixth method, one of the collinear x^0 x^n x^np1 distributions.
  // First choose n.
  temp = r.random(); 
  for(n = 1; double(n)/double(nverts) < temp; n++){} 
  // This exits with n = first n with n/nverts >= temp.
  np1 = cyclic(n+1,nverts);
  // Now choose xs.
  xs = xsmax*pow(r.random(),1.0/cparam);
  // Now choose x_i for i not in S.
	for(i = 0; i < nverts - 3; i++) rval[i] = r.random(); // Get nverts - 3 random numbers in (0,1).
	DESsort(rval,nverts - 3);                             // Sort the random numbers.
  rval[nverts - 3] = 1.0;                               // An artificial largest r.
  j = -1;
	for(i = 1; i <= nverts; i++)
	  {
      if(i != n && i != np1)
        {
        j = j + 1;
        if (j == 0)
          {
          x[i] = (1.0 - xs)*rval[j];
          }
        else
          {
          x[i] = (1.0 - xs)*(rval[j] - rval[j-1]);
          }
        }
	  }
  // Now choose x_i for i in S.
  rval[0] = r.random(); // Get 2 random numbers in (0,1).
  rval[1] = r.random();
  DESsort(rval,2);      // Sort the random numbers.
  x[0] = xs*rval[0];
  x[n] = xs*(rval[1] - rval[0]);
  x[np1] = xs*(1.0 - rval[1]);
  OK = true;
  }
else if(method == 7)
  {
  // Seventh method, one of the collinear-soft distributions.
  // First choose n.
  temp = r.random(); 
  for(n = 1; double(n)/double(nverts) < temp; n++){} 
  // This exits with n = first n with n/nverts >= temp.
  nm1 = cyclic(n-1,nverts);
  np1 = cyclic(n+1,nverts);
  np2 = cyclic(n+2,nverts);
  // We generate the needed information for either collinear-soft(-) or collinear-soft(+).
  // Notation follows the collinear-soft(-) case.
  // Find y(0)
  y0 = r.random();
  // Find y(j) for j != 0, n-1, n, n+1 (or 0,n,n+1,n+2), 
  // which we label as yy[i] for i = 0,...,nverts-4.
  if(nverts < 5)
	{
	if(nverts == 4) yy[0] = 1.0;
	else assert(nverts >= 4) ;
	}
  else // nverts is at least 5;
	{
	for(i = 0; i < nverts - 4; i++) rval[i] = r.random(); // Get nverts - 4 random numbers in (0,1).
	DESsort(rval,nverts - 4);                              // Sort the random numbers.
	yy[0] = rval[0];
	for(i = 1; i <= nverts - 5; i++)
	  {
	  yy[i] = rval[i] - rval[i-1];
	  }
	yy[nverts - 4] = 1.0 - rval[nverts-5];
	}
  // Now, will it be collinear-soft(+) or collinear-soft(-)?
  temp = r.random(); 
  if(temp > 0.5) // We look for collinear-soft(-) choice.
    {
    // Start with xA, xB, and x[n].
    xA = xAminusmax[n]*pow(r.random(),1.0/cparam);
    xB = xBminusmax[n]*sqrt(r.random());
    x[n] = xnminusmax[n]*r.random();
    x[0] = 2.0*x[n]*xA*xB*y0;
    x[np1] = x[n]*xA*xB*(1.0 - y0);
    // Find x[cyclic(n+2,nverts),..., x[cyclic(n-2,nverts)].
    for(i = 0; i <= nverts - 4; i++)
      {
      temp = abs(s[nm1][np1]/s[cyclic(n+2+i,nverts)][n]);
      x[cyclic(n+2+i,nverts)] = xA * temp * yy[i];
      }
    // Finally, find x[n-1].
    temp = 0.0;
    for(i = 0; i <= nverts; i++)
      {
      if(i != nm1) temp = temp + x[i];
      }
      OK = false;
    if(temp < 1.0) OK = true;
    x[nm1] = 1.0 - temp;
  	}
  else  // We look for collinear-soft(+) choice.
    {
    // Start with xA, xB, and x[n+1].
    xA = xAplusmax[n]*pow(r.random(),1.0/cparam);
    xB = xBplusmax[n]*sqrt(r.random());
    x[np1] = xnp1plusmax[n]*r.random();
    x[0] = 2.0*x[np1]*xA*xB*y0;
    x[n] = x[np1]*xA*xB*(1.0 - y0);
    // Find x[cyclic(n+3,nverts),..., x[cyclic(n-1,nverts)].
    for(i = 0; i <= nverts - 4; i++)
      {
      temp = abs(s[n][np2]/s[cyclic(n+3+i,nverts)][np1]);
      x[cyclic(n+3+i,nverts)] = xA * temp * yy[i];
      }
    // Finally, find x[n+2].
    temp = 0.0;
    for(i = 0; i <= nverts; i++)
      {
      if(i != np2) temp = temp + x[i];
      }
      OK = false;
      if(temp < 1.0) OK = true;
      x[np2] = 1.0 - temp;
    } // Close if(temp > 0.5) ... else ...
  }
else if(method == 8)
  {
  // Eigth method, the double parton scattering distribution.
  if(dps.n == 0)
    {
    cout << "This cannot happen. We got dps.n = 0 when it must be 1." << endl;
    exit(1);
    }
  else // dps.n = 1 so we should use this method to choose a point.
    {
    // First, the variables y^j for j in S, chosen uniformly so they sum to 1.
    for(i = 0; i < nverts - 4; i++) rval[i] = r.random(); // Get nverts - 4 random numbers in (0,1).
    DESsort(rval,nverts - 4);                             // Sort the random numbers.
    j = 0;
    for(i = 1; i <= nverts; i++)
      {
      if(i != dps.A && i != dps.Ap1 && i != dps.B && i != dps.Bp1)
        {
        j = j + 1;
        if(j == 1) yy[i] = rval[0];
        else yy[i] = rval[j - 1] - rval[j - 2];
        }
      }
    yy[0] = 1.0 - rval[j - 1];
    }
    // Now, the other integration variables.
    rbar = r.random();
    t = r.random();
    u = r.random();
    tildeomega = r.random();
    if(r.random() < 0.5) signA = 1.0;
    else signA = -1.0;
    if(r.random() < 0.5) signB = 1.0;
    else signB = -1.0;
    // From these, we generate other variables.
    axbar = 0.5;                                           // Must match choice in DESxpoint::rho.
    xbarpower = 1.0/(1.0 - axbar);
    if(rbar < 0.5)
      {
      xbar = 0.5*pow(2.0*rbar,xbarpower);
      oneminusxbar = 1.0 - xbar;
      }
    else
      {
      oneminusxbar = 0.5*pow(2.0*(1.0 - rbar),xbarpower);
      xbar = 1.0 - oneminusxbar;
      }
    deltau = 0.1*xbar*oneminusxbar*abs(dps.df);             // Must match choice in DESxpoint::rho.
    omegafactorA = 1.0 + signA*(2.0*dps.fA - 1.0);
    omegafactorB = 1.0 + signB*(2.0*dps.fB - 1.0);
    xumax = 1.0*xbar*oneminusxbar*omegafactorA*omegafactorB; // Must match choice in DESxpoint::rho.
    xu = 1.0 + pow(xumax/deltau, nverts - 2);
    xu = pow(xu,t);
    xu = xu - 1.0;
    xu = pow(xu, 1.0/(nverts - 2.0));
    xu = deltau*xu;
    temp = pow(u,1.0/(nverts - 3.0));
    xs = xu*temp;
    xt = xu*(1.0 - temp);
    oneminusxs = 1.0 - xs;
    sqrtxt = sqrt(xt);
    omegamax =   log(0.5*omegafactorA*xbar/sqrtxt);          // Must match choice in DESxpoint::rho.
    omegamin = - log(0.5*omegafactorB*(1.0 - xbar)/sqrtxt);  // Must match choice in DESxpoint::rho.
    omega = omegamin + (omegamax - omegamin)*tildeomega;
    expomega = exp(omega);
    yA = dps.fA*xbar - signA*sqrtxt*expomega;
    yAp1 = dps.oneminusfA*xbar + signA*sqrtxt*expomega;
    yB = dps.fB*oneminusxbar - signB*sqrtxt/expomega;
    yBp1 = dps.oneminusfB*oneminusxbar + signB*sqrtxt/expomega;
    // Now we construct the x^i.
    for(i = 0; i <= nverts; i++)
      {
      if(i != dps.A && i != dps.Ap1 && i != dps.B && i != dps.Bp1)
        {
        x[i] = xs*yy[i];
        }
      }
    x[dps.A]   = oneminusxs*yA;
    x[dps.Ap1] = oneminusxs*yAp1;
    x[dps.B]   = oneminusxs*yB;
    x[dps.Bp1] = oneminusxs*yBp1;
    // Check the point.
    OK = true;
    sum = 0.0;
    for(i = 0; i <= nverts; i++)
      {
      if(x[i] < 0.0) OK = false;
      sum = sum + x[i];
      }
    if(abs(sum - 1.0) > tiny)
      {
      cout << "Sum is wrong for method 8 "<< sum << endl;
      for(i = 0; i <= nverts; i++) cout << "x["<<i<<"] = "<< x[i] << endl;
      exit(1);
      }
  }
else if(method == 9)
  {
  // Ninth method, one of the collinear distribtions.
  if(collN == 0)
    {
    cout << "This cannot happen. We got collN = 0 when it must be > 0." << endl;
    exit(1);
    }
  else // collN > 0 so we should use this method to choose a point.
    {
    // First choose the index n corresponding to which singularity we want.
    temp = r.random(); 
    for(n = 1; double(n)/double(collN) < temp; n++){} 
    // This exits with n = first n with n/collN >= temp.
    // First, the variables y^j for j in S, chosen uniformly so they sum to 1.
    for(i = 0; i < nverts - 4; i++) rval[i] = r.random(); // Get nverts - 4 random numbers in (0,1).
    DESsort(rval,nverts - 4);                             // Sort the random numbers.
    j = 0;
    for(i = 1; i <= nverts; i++)
      {
      if(i != coll[n].N && i != coll[n].Np1 && i != coll[n].A && i != coll[n].B)
        {
        j = j + 1;
        if(j == 1) yy[i] = rval[0];
        else yy[i] = rval[j - 1] - rval[j - 2];
        }
      }
    yy[0] = 1.0 - rval[j - 1];
    // Now, the other integration variables.
    tu     = r.random();
    xa     = r.random();
    u      = r.random();
    rcoll  = r.random();
    // Now we construct the x^i.
    lambdaR = 1/sqrt(coll[n].df); // This choice matches the choice in DESxpoint::rho.
    acoll = 1.0*coll[n].df;       // This choice matches the choice in DESxpoint::rho. 
    xu = tu*acoll/(1.0 + acoll - tu);
    temp = pow(u,1.0/(nverts - 3.0));
    xs = xu*temp;
    oneminusxs = 1.0 - xs;
    ycoll = xu*(1.0 - temp);
    dfsq = pow(coll[n].df,2);
    xt = pow(sqrt(ycoll)*rcoll/lambdaR/coll[n].df,twothirds);
    oneminusxt = 1.0 - xt;
    xbarsq = pow(ycoll*lambdaR*coll[n].df/rcoll,twothirds) - dfsq;
    if(xbarsq <= 0.0 || oneminusxt <= 0.0)
      {
      for(i = 0; i <= nverts; i++) xOUT.in(i) = 0.0; 
      xOK = false;
      return;       // A bad point.
      }
    if(r.random() < 0.5) xbar = sqrt(xbarsq);
    else xbar = - sqrt(xbarsq);
    yN = oneminusxt*(coll[n].f - xbar);
    yNp1 = oneminusxt*(1.0 - coll[n].f + xbar);
    yA = xt*xa;
    yB = xt*(1.0 - xa);
    for(i = 0; i <= nverts; i++)
      {
      if(i != coll[n].N && i != coll[n].Np1 && i != coll[n].A && i != coll[n].B)
        {
        x[i] = xs*yy[i];
        }
      }
    x[coll[n].N]   = oneminusxs*yN;
    x[coll[n].Np1] = oneminusxs*yNp1;
    x[coll[n].A]   = oneminusxs*yA;
    x[coll[n].B]   = oneminusxs*yB;
    // Check the point.
    OK = true;
    sum = 0.0;
    for(i = 0; i <= nverts; i++)
      {
      if(x[i] < 0.0) OK = false;
      sum = sum + x[i];
      }
    if(abs(sum - 1.0) > tiny)
      {
      cout << "Sum is wrong for method 9 "<< sum << endl;
      for(i = 0; i <= nverts; i++) cout << "x["<<i<<"] = "<< x[i] << endl;
      exit(1);
      }
    }
  }
else
  {
  cout << "Non-existent method for choosing point detected. "<< method << endl;
  exit(1);
  }
// End various methods for choosing the point.
//------
// Check that we are not too near a boundary.
// Set xOUT = x and xOK = OK.
for(i = 0; i <= nverts; i++)
  {
  if (x[i] < xcutoff) OK = false;
  xOUT.in(i) = x[i]; 
  }
xOK = OK;
}
//--------------------------------------------------------------------------------------------------
// Return rho(x) for a point x.
double DESxpoint::rho()
{
static double rhoval;
static int n,i,nm1,np1,np2;
static double gval,capGval,xbar,omegamax,omegamin;
static double partialrho;
static double probability;
static double xs,rtxs,svalue,ys,yn1;
static double xA,xB;
static double temp,prefactor,xAfactor;
static const double tiny = 1.0e-40;
static double denom1,denom2;
static double deltau,xumax,oneminusxbar;
static double xt,xu,acoll;
static double axbar;
static double r,lambdaR;
static double signA,signB,omegafactorA,omegafactorB;
static double tempA,tempB;
static double sqrtxt;
// Contribution from the uniform distribution, method 1.
probability = alpha[1];
rhoval = probability*double(factorial(nverts));
// Contribution distribution with small x^0, uniform x^i distribution, method 10.
probability = alpha[10];
temp = pow(1.0 - x[0],nverts);
partialrho = probability*double(factorial(nverts))*c10/pow(1.0 - temp,1.0 - c10);
rhoval = rhoval + partialrho;
// Contribution distribution with large x^0, uniform x^i distribution, method 11.
probability = alpha[11];
if(probability > tiny)
  {
  partialrho = probability*double(factorial(nverts-1))*c11/pow(1.0 - x[0],nverts - c11);
  rhoval = rhoval + partialrho;
  }
// The contributions for methods 2,..., 7 depend on index n, so we loop over n.
for(n = 1; n <= nverts; n++)
  {
  nm1 = cyclic(n-1,nverts);
  np1 = cyclic(n+1,nverts);
  np2 = cyclic(n+2,nverts);
  // Contributions from soft-n distributions, method 2.
  probability = alpha[2];
  // For instruction if (method == 8 && dps.n == 0) method = 2; 
  // Realocate to soft-n method if dpf method not needed.
  if (dps.n == 0) probability = probability + alpha[8];
  // For instruction if (method == 9 && collN == 0) method = 2; 
  // Realocate to soft-n method if coll method not needed.
  if (collN == 0) probability = probability + alpha[9];
  //
  gval = g(n);
  if(gval < capA[n])
    {
    xbar = x[nm1]*x[np1];
    omegamax = log(1/sqrt(xbar));
    capGval = gval/lambda[n];
    partialrho = probability/wsoftsum*cparam*double(factorial(nverts - 2))/2.0;
    partialrho = partialrho/omegamax/pow(capGval,nverts - 1.0 - cparam);
    rhoval = rhoval + partialrho;
    }
  // Contributions from modified collinear distributions, method 3.
  probability = alpha[3];
  if(probability > tiny)
    {
    xbar = x[0] + x[n] + x[np1];
    partialrho = probability/double(nverts)*double(factorial(nverts - 2))*2.0/pow(xbar,2);
    rhoval = rhoval + partialrho;
    }
  // Contributions from modified soft distributions, method 4.
  probability = alpha[4];
  if(probability > tiny)
    {
    xbar = x[nm1] + x[n] + x[np1];
    partialrho = probability/double(nverts)*double(factorial(nverts - 2))*2.0/pow(xbar,2);
    rhoval = rhoval + partialrho;
    }
  // Contributions from the "x0" collinear distritutions, method 5.
  probability = alpha[5];
  if(probability > tiny)
    {
    // First, contribution for x[n-1] & x[n] large.
    if(s[nm1][np1] > 0)
      {
      xs = 0;
      for(i = n + 1; i <= n + nverts - 2; i++)
        {
        xs = xs + x[cyclic(i,nverts)];
        }
      rtxs = sqrt(xs);
      svalue = s[nm1][np1];
      ys = svalue/(svalue  + 2.0*m0sq);
      omegamin = atan(- ys/rtxs);
      omegamax = atan((1.0 - ys)/rtxs);
      denom1 = 1.0 - x[0] - xs;
      denom2 = 1.0 - xs;
      if(denom1 < tiny) denom1 = tiny;
      if(denom2 < tiny) denom2 = tiny;
      yn1 = x[n]/denom1;
      partialrho = probability/Ncollx0;
      partialrho = partialrho*cparam*double(factorial(nverts - 3))*pow(xs,cparam + 2 - nverts);
      partialrho = partialrho*rtxs*(1.0 - xs + rtxs)/denom2/pow((x[0] + rtxs),2);
      partialrho = partialrho*rtxs/(pow(yn1 - ys,2) + xs);
      partialrho = partialrho/denom1/(omegamax - omegamin);
      rhoval = rhoval + partialrho;
      }
    // Second, contrbution for x[n+1] & x[n+2] large.
    if(s[n][np2] > 0)
      {
      xs = 0;
      for(i = n + 3; i <= n + nverts; i++)
        {
        xs = xs + x[cyclic(i,nverts)];
        }
      rtxs = sqrt(xs);
      svalue = s[n][np2];
      ys = svalue/(svalue + 2.0*m0sq);
      omegamin = atan(- ys/rtxs);
      omegamax = atan((1.0 - ys)/rtxs);
      denom1 = 1.0 - x[0] - xs;
      denom2 = 1.0 - xs;
      if(denom1 < tiny) denom1 = tiny;
      if(denom2 < tiny) denom2 = tiny;
      yn1 = x[np1]/denom1;
      partialrho = probability/Ncollx0;
      partialrho = partialrho*cparam*double(factorial(nverts-3))*pow(xs,cparam + 2 - nverts);
      partialrho = partialrho*rtxs*(1.0 - xs + rtxs)/denom2/pow((x[0] + rtxs),2);
      partialrho = partialrho*rtxs/(pow(yn1 - ys,2) + xs);
      partialrho = partialrho/denom1/(omegamax - omegamin);
      rhoval = rhoval + partialrho;
      }
    } // Close if(probability > tiny).
  // Contributions form the collinear x^0 x^n x^np1 distributions, method 6.
  probability = alpha[6];
  if(probability > tiny)
    {
    xs = x[0] + x[n] + x[np1];
    if(xs < xsmax)
      {
      partialrho = probability/double(nverts)
                  *2.0*double(factorial(nverts - 3))*cparam/pow(xsmax,cparam)
                  *pow(xs,cparam - 3.0)*pow(1.0 - xs,3 - nverts);
      rhoval = rhoval + partialrho;
      }
    }
  // Contributions from the collinear-soft(+/-) distributions, method 7.
  probability = alpha[7];
  if(probability > tiny)
    {
    prefactor = probability/double(2*nverts)*cparam*double(factorial(nverts - 4));
    // First, contributions from the collinear-soft(-) distributions.
    xA = 0.0;
    xAfactor = 1.0;
    for(i = 1; i <= nverts; i++)
      {
    if(i != nm1 && i != n && i != np1)
      {
      temp = abs(s[i][n]/s[nm1][np1]);
      xA = xA + x[i] * temp;
      xAfactor = xAfactor * temp;
      }
    }
    xB = (x[np1] + 0.5*x[0])/x[n]/xA;
    if(xA < xAminusmax[n] && xB < xBminusmax[n] && x[n] < xnminusmax[n])
      {
      // xA, xB, and x[n] are in range, so we have a nonzero partialrho.
      partialrho = prefactor
                /pow(xAminusmax[n],cparam)/pow(xBminusmax[n],2)/xnminusmax[n]
                *xAfactor;
      partialrho = partialrho /pow(x[n],2) /pow(xA,nverts - cparam - 1.0);
      rhoval = rhoval + partialrho;
      }
    // Next, contributions from the collinear-soft(+) distributions.
    // prefactor = pcsoft/double(2*nverts)*cparam*factorial(nverts - 4) as above
    xA = 0.0;
    xAfactor = 1.0;
    for(i = 1; i <= nverts; i++)
      {
      if(i != np2 && i != n && i != np1)
      {
      temp = abs(s[i][np1]/s[np2][n]);
      xA = xA + x[i] * temp;
      xAfactor = xAfactor * temp;
      }
      }
    xB = (x[n] + 0.5*x[0])/x[np1]/xA;
    if(xA < xAplusmax[n] && xB < xBplusmax[n] && x[np1] < xnp1plusmax[n])
      {
      // xA, xB, and x[n+1] are in range, so we have a nonzero partialrho.
      partialrho = prefactor
                /pow(xAplusmax[n],cparam)/pow(xBplusmax[n],2)/xnp1plusmax[n]
          *xAfactor;
      partialrho = partialrho /pow(x[np1],2) /pow(xA,nverts - cparam - 1.0);
      rhoval = rhoval + partialrho;
      }
    } // Close if(probability > tiny).
  } // Close for(n = 1; n <= nverts; n++).
// Now there are two more for which we do not loop over n.
//
// Contribution from the double parton scattering distribution, method 8.
probability = alpha[8];
if(probability > tiny && dps.n > 0)
  {
  axbar = 0.5;                                              // Must match choice in DESxpoint::neu.
  prefactor = probability*0.25*(1.0 - axbar)*double(factorial(nverts - 2));
  xs = x[0];
  for(i = 1; i <= nverts; i++)
    {
    if(i != dps.A && i != dps.Ap1 && i != dps.B && i != dps.Bp1) xs = xs + x[i];
    }
  xbar = (x[dps.A] + x[dps.Ap1])/(1.0 - xs);
  oneminusxbar = 1.0 - xbar;
  tempA = dps.fA*x[dps.Ap1] - dps.oneminusfA*x[dps.A];
  tempB = dps.fB*x[dps.Bp1] - dps.oneminusfB*x[dps.B];
  if (tempA > 0) signA = 1.0;
  else signA = - 1.0;
  if (tempB > 0) signB = 1.0;
  else signB = - 1.0;
  xt = abs(tempA*tempB) /pow(1.0 - xs,2);
  sqrtxt = sqrt(xt);
  xu = xt + xs;
  deltau = 0.1*xbar*oneminusxbar*abs(dps.df);               // Must match choice in DESxpoint::neu.
  omegafactorA = 1.0 + signA*(2.0*dps.fA - 1.0);
  omegafactorB = 1.0 + signB*(2.0*dps.fB - 1.0);
  xumax = 1.0*xbar*oneminusxbar*omegafactorA*omegafactorB;  // Must match choice in DESxpoint::neu.
  omegamax =   log(0.5*omegafactorA*xbar/sqrtxt);           // Must match choice in DESxpoint::neu.
  omegamin = - log(0.5*omegafactorB*(1.0 - xbar)/sqrtxt);   // Must match choice in DESxpoint::neu.
  if (xu < xumax)
    {
    partialrho = prefactor/pow(2.0*min(xbar,oneminusxbar),axbar);
    partialrho = partialrho/pow(1.0 - xs,3)/(omegamax - omegamin);
    partialrho = partialrho/log(1.0 + pow(xumax/deltau,nverts-2));
    partialrho = partialrho/(pow(xu,nverts-2) + pow(deltau,nverts-2));
    rhoval = rhoval + partialrho;
    }
  }
// Contributions from the collinear distribtions, method 9.
probability = alpha[9];
if(probability > tiny && collN > 0)
  {
  prefactor = probability*1.5*double(factorial(nverts - 3))/double(collN);
  for(n = 1; n <= collN; n++)
    {
    lambdaR = 1/sqrt(coll[n].df);  // This choice matches with DESxpoint::neu.
    acoll = 1.0*coll[n].df;        // This choice matches the choice in DESxpoint::neu.
    xs = x[0];
    for(i = 1; i <= nverts; i++)
      {
      if(i != coll[n].N && i != coll[n].Np1 && i != coll[n].A && i != coll[n].B) xs = xs +x[i];
      }
    xt = (x[coll[n].A] + x[coll[n].B])/(1.0 - xs);
    xbar = (coll[n].f*x[coll[n].Np1] - (1.0 - coll[n].f)*x[coll[n].N])/(1.0 - xs)/(1.0 - xt);
    temp = pow(xbar,2) + pow(coll[n].df,2);
    xu = temp*xt + xs;
    r = lambdaR*coll[n].df*xt/sqrt(temp);
    if(xu < 1.0 && r < 1.0)
      {
      partialrho = prefactor/(1.0 - xt)/pow(1.0 - xs,3);
      partialrho = partialrho*abs(xbar)*lambdaR*coll[n].df/sqrt(temp);
      partialrho = partialrho*acoll*(1.0 + acoll)/pow(xu + acoll,2);
      partialrho = partialrho/pow(xu,nverts - 3);
      rhoval = rhoval + partialrho;
      }
    }
  }
//
return rhoval;
}
//--------------------------------------------------------------------------------------------------
// Reset x and get rho(x) from rho(). 
double DESxpoint::altrho(const DESxarray& xIN)
{
static double rhoval;
for(int i = 0; i <= nverts; i++)
  {
  x[i] = xIN(i);
  }
rhoval = rho();
return rhoval;
}
//--------------------------------------------------------------------------------------------------
// Scaled version of approximation to \Lambda in soft region n
// g_n(x). Note that x, a_{ni} and b_n are known to class DESxpoint.
//
double  DESxpoint::g(const int& n)
{
static int i,ii;
static double gg;
gg = b[n]*x[cyclic(n-1,nverts)]*x[cyclic(n+1,nverts)];
for(i = 2; i <=  nverts - 2; i++)
  {
  ii = cyclic(n+i,nverts);
  gg = gg + a[n][ii]*x[ii];
  }
gg = gg + a[n][0]*x[0];
return gg;
}
// End of implementation code for DESxpoint.
//==================================================================================================
// Implementation code for class DESxdeform.
//
// Constructor
DESxdeform::DESxdeform(int nvertsIN, double rtssqIN, double lambdadeformIN)
{
nverts = nvertsIN;
rtssq = rtssqIN;
lambdadeform = lambdadeformIN;
for(int i = 0; i <= nverts; i++)
  {
  for(int j = 0; j <= nverts; j++) s[i][j] = 0.0;
  }
}
//--------------------------------------------------------------------------------------------------
// Sets private Smatrix and other private variables.
void DESxdeform::newS(DES_Smatrix sIN)
{
static int i,j;
for(i = 0; i <= nverts; i++)
  {
  for(j = 0; j <= nverts; j++) s[i][j] = sIN(i,j);
  }
assert(nverts == sIN.getnverts());
assert(rtssq == sIN.getrtssq());
}
//--------------------------------------------------------------------------------------------------
// Returns deformed point z = x + i y and deformation jacobian dz/dx.
void DESxdeform::deform(/*in*/ DESxarray& x, 
                       /*out*/ DESzarray& z, complex<double>& jacobian, double& etasum)
{
static const complex<double>ii(0.0,1.0);
static int i,j,nn,nnp1,nnp2,nnm1;
static double yy[maxsize];                      // Deformation varialbe eta[i].
static double w[maxsize];                       // Gradient of Lambda.
static complex<double> dzdx[maxsize][maxsize];  // Matrix dz(i+1)/dx(j+1) for i,j in {1,...,nverts}.
static double dydx[maxsize][maxsize];           // Matrix dy(i)/dx(j) for i,j in {0,...,nverts}.
static double deformparam,c_g;
static double temp;
static double dsumdxj;
static double g,h,gplus,gminus,dgplusdxj,dgminusdxj;
static double dgdx[maxsize],dloghdx[maxsize];
static double deltaS;
static int iplus,iminus;
//
deformparam = 6.0*lambdadeform/rtssq;  // For standard method.  Also used for extra deformation.
c_g = 0.5;                             // For extra deformation in collinear direction.
//
//=================
// The standard deformation method. Applies to all terms.
// This defines w[i] and gives the starting values for yy[i] and dydx[i][j].
//
w[0] = 0.0;                     // We don't actually use this.
for(i = 1; i <= nverts; i++)
  {
  w[i] = 0.0;
  for(j = 1; j <= nverts; j++)
    {
    w[i] = w[i] + s[i][j]*x(j);
    }
  }
yy[0] = 0.0;                     // We define this to vanish.
for(i = 1; i <= nverts; i++)
  {
  yy[i] = deformparam*x(i)*w[i];
  }
// Construct the matrix dydx.
for(j = 0; j  <= nverts; j++)
  {
  dydx[0][j] = 0.0;              // We don't actually use this.
  }
for(i = 1; i <= nverts; i++)
  {
  dydx[i][0] = 0.0;              // We don't actually use this.
  for(j = 1; j  <= nverts; j++)
    {
    dydx[i][j] = deformparam*x(i)*s[i][j];
    }
    dydx[i][i] = dydx[i][i] + deformparam*w[i];
  }
//==== Extra deformation along direction of collinear singularity =======
for(nn = 1; nn <= nverts; nn++)
  {
  nnp1 = nn + 1;
  if(nnp1 > nverts) nnp1 = nnp1 - nverts;
  nnp2 = nn + 2;
  if(nnp2 > nverts) nnp2 = nnp2 - nverts;
  nnm1 = nn - 1;
  if(nnm1 < 1) nnm1 = nnm1 + nverts;
  // Check that point x is OK for this nn.
  if(x(nn) + x(nnp1) > 0.5 && x(nn) > x(nnp2) && x(nnp1) > x(nnm1))
    {
    h = (2.0*x(nn) + 2.0*x(nnp1) - 1.0)*4.0*(x(nn) - x(nnp2))*(x(nnp1) - x(nnm1));
    // Derivatives of log(h(x)).
    for(j = 0; j <= nverts; j++) dloghdx[j] = 0.0;  // We don't actually use dloghdx[0].
    dloghdx[nn]   =  2.0/(2.0*x(nn) + 2.0*x(nnp1) - 1.0) + 1.0/(x(nn) - x(nnp2));
    dloghdx[nnp1] =  2.0/(2.0*x(nn) + 2.0*x(nnp1) - 1.0) + 1.0/(x(nnp1) - x(nnm1));
    dloghdx[nnp2] = - 1.0/(x(nn) - x(nnp2));
    dloghdx[nnm1] = - 1.0/(x(nnp1) - x(nnm1));
    //
    gplus = lambdadeform/deformparam;
    gminus = lambdadeform/deformparam;
    iplus = 0;
    iminus = 0;
    for(i = 1; i <= nverts; i++)
      {
      if ( i != nn && i != nnp1 )  // i is in S.
        {
        deltaS = s[i][nnp1] - s[i][nn];
        if(deltaS > 0) // i is in S+.
          {
          temp = pow(w[i],2)/deltaS;
          if(temp < gplus)
            {
            gplus = temp;
            iplus = i;
            }
          }
        if(deltaS < 0) // i is in S-.
          {
          temp = - pow(w[i],2)/deltaS;
          if(temp < gminus)
            {
            gminus = temp;
            iminus = i;
            }
          }
        } // Close if ( i != nn && i != nnp1 ).
      } // Close for( i = 1; i <=n; i++).
    g = h*c_g*(gplus  - gminus);          // Parameter c_g was set at top. Default 0.5.
    // Calculate derivatives.
    dgdx[0] = 0.0;                // We don't actually use this.
    if(iplus != 0 && iminus != 0) // So g+ and g- are both non-trivial.
      {  
      for(j = 1; j <= nverts; j++)
        {
        dgplusdxj  = 2.0*w[iplus]/(s[iplus][nnp1] - s[iplus][nn])*s[iplus][j];
        dgminusdxj = 2.0*w[iminus]/(s[iminus][nn] - s[iminus][nnp1])*s[iminus][j];
        dgdx[j] = g * dloghdx[j] + c_g*h*(dgplusdxj - dgminusdxj);
        }
      }
    else if(iplus != 0) // So g- is just a constant.
      {
      for(j = 1; j <= nverts; j++)
        {
        dgplusdxj  = 2.0*w[iplus]/(s[iplus][nnp1] - s[iplus][nn])*s[iplus][j];
        dgminusdxj = 0.0;
        dgdx[j] = g * dloghdx[j] + c_g*h*(dgplusdxj - dgminusdxj);
        }
      }
    else if(iminus != 0) // So g+ is just a constant.
      {
      for(j = 1; j <= nverts; j++)
        {
        dgplusdxj  = 0.0;
        dgminusdxj = 2.0*w[iminus]/(s[iminus][nn] - s[iminus][nnp1])*s[iminus][j];
        dgdx[j] = g * dloghdx[j] + c_g*h*(dgplusdxj - dgminusdxj);
        }
      }
    else // Both g+ and g- are consants.
      {
      for(j = 1; j <= nverts; j++)
        {
        dgminusdxj = 0.0;
        dgplusdxj  = 0.0;
        dgdx[j] = g * dloghdx[j] + c_g*h*(dgplusdxj - dgminusdxj);
        }
      }
    // We now have g and dloggdx. Calculate deformation yy and its derivatives.
    yy[nn] = yy[nn] + deformparam * g;
    yy[nnp1] = yy[nnp1] - deformparam * g;
    for(j = 1; j <= nverts; j++)
      {
      dydx[nn][j]   = dydx[nn][j]   + deformparam * dgdx[j];
      dydx[nnp1][j] = dydx[nnp1][j] - deformparam * dgdx[j];
      }
    } // Close if(x(nn) + x(nnp1) > 0.5 && x(nn) > x(nnp2) && x(nnp1) > x(nnm1))
  } // Close for(nn = 1; nn <= nverts; nn++)
//
//-------------
// We have now specified the deformation yy. 
// We set z = (x + i yy)/(1 + sum yy)
//
etasum = 0.0;
for(i = 1; i <= nverts; i++) etasum = etasum + yy[i];
for(i = 0; i <= nverts; i++) // Note that this includes definition of z(0).
  {
  z.in(i) = (x(i) + ii*yy[i])/(1.0 + ii*etasum);
  }
// The matrix dz/dx.
// In principle, here we treat x(0) as a dependent variable.
// However, since there is no x(0) dependence, that does not matter.
// The independent variables are x(1),...,x(nverts).
// Our labeling is that what any reasonable person
// would call dzdy(i,j) we call dzdy[i-1][j-1].
for(j = 1; j <= nverts; j++)
  {
  dsumdxj = 0.0;
  for(i = 1; i <= nverts; i++)
    {
    dsumdxj = dsumdxj + dydx[i][j];
    }
  for(i = 1; i <= nverts; i++)
    {
    dzdx[i-1][j-1] = ii*( dydx[i][j] - z(i)*dsumdxj)
                     /(1.0 + ii*etasum);
    }
  dzdx[j-1][j-1] =  dzdx[j-1][j-1] + 1.0/(1.0 + ii*etasum);
  }
// Finally, we calculate the determinant of this matrix.
jacobian = Determinant(dzdx,nverts);
}
// End of implementation code for DESxdeform.
//==================================================================================================
// Implementation code for class DESellpoint.
//
// Constructor.
DESellpoint::DESellpoint(int seedIN, int nvertsIN)
  : r(seedIN)     // Construct random number generator instance r.
{
seed = seedIN;
nverts = nvertsIN;
Aradial = double(nverts)/2.0 - 0.5;  // Presumably this is a good choice.
}
// There is no default constructor.
//--------------------------------------------------------------------------------------------------
// Generate new point ell.
// Here ellOUT is a (complex) vector {ell[0], ell[1], ell[2], ell[3]}.
void DESellpoint::neu(DEScmplx4vector& ellOUT)
{
// First construct point at random on unit sphere in 4 dimensions.
// d^3 \Omega = (1/2) d\theta_A d\theta_B d cos^2(theta_C) where
// 0 < theta_A < 2 Pi, 0 < theta_A < 2 Pi, 0 < theta_C < Pi/2 and
// ell[0] = cosC*cosA; ell[1] = cosC*sinA; ell[2] = sinC*cosB; ell[3] = sinC*sinB;
// The density of points is 1/area = 1/(2 Pi^2).
static double cosA,sinA,cosB,sinB,cosC,sinC,thetaA,thetaB;
static double temp;
static int mu;
static const double twoPI = 6.283185307179586;
thetaA = r.random()*twoPI;
sinA = sin(thetaA);
cosA = cos(thetaA);
thetaB = r.random()*twoPI;
sinB = sin(thetaB);
cosB = cos(thetaB);
temp = r.random();
cosC = sqrt(temp);
sinC = sqrt(1.0 - temp);
ell.in(0) = cosC*cosA;
ell.in(1) = cosC*sinA;
ell.in(2) = sinC*cosB;
ell.in(3) = sinC*sinB;
//
// Now fix the length of ell.
ellEsq = pow(r.random(),-1.0/(Aradial - 1.0)) - 1.0;
temp = sqrt(ellEsq);
for(mu = 0; mu <= 3; mu++)
  {
   ell.in(mu) = temp*ell(mu);
   ellOUT.in(mu) = ell(mu); // This converts to complex.
   }
}
//--------------------------------------------------------------------------------------------------
// Return rho(ell) for a point ell.
double DESellpoint::rho()
{
static double rhoval;
static const double invtwoPIsq = 0.05066059182116889;  // 1/(2*Pi^2)
rhoval = invtwoPIsq;
rhoval = rhoval * 2.0*(Aradial - 1.0)/ellEsq/pow((1.0 + ellEsq),Aradial);
return rhoval;
}
// End of implementation code for DESellpoint.
//==================================================================================================
// Implementation code for class DESreal4vector.

// Default constructor
DESreal4vector::DESreal4vector()
{
for(int mu = 0; mu <= 3; mu++)
  {
  vector[mu] = 0.0;
  }
}
//--------------------------------------------------------------------------------------------------
// v.in(mu) for input, v.in(mu) = r
double&  DESreal4vector::in(const int& mu)
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// v(mu) for output,  r = v(mu)
const double&  DESreal4vector::operator()(const int& mu) const
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// Sum of two vectors: v3 = v1 + v2
DESreal4vector DESreal4vector::operator+(const DESreal4vector& v2) const
{
DESreal4vector sum;
for(int mu = 0; mu <= 3; mu++)
  {
  sum.vector[mu] =  vector[mu] + v2.vector[mu];
  }
return sum;
}
//--------------------------------------------------------------------------------------------------
// Difference of two vectors: v3 = v1 - v2
DESreal4vector DESreal4vector::operator-(const DESreal4vector& v2) const
{
DESreal4vector difference;
for(int mu = 0; mu <= 3; mu++)
  {
  difference.vector[mu] =  vector[mu] - v2.vector[mu];
  }
return difference;
}
//--------------------------------------------------------------------------------------------------
// Vector times scalar: v2 = v1*c
DESreal4vector  DESreal4vector::operator*(const double& a) const
{
DESreal4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.vector[mu] = a * vector[mu];
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// -1 times vector: v2 = - v1
DESreal4vector  DESreal4vector::operator-() const
{
DESreal4vector negative;
for(int mu = 0; mu <= 3; mu++)
  {
  negative.vector[mu] = - vector[mu];
  }
return negative;
}

// Scalar times vector: v2 = c*v1
DESreal4vector operator*(const double& a, const DESreal4vector& v)
{
DESreal4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.in(mu) = a * v(mu);  // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Dot product of two 4-vectors: dot(v1,v2)
double dot(const DESreal4vector& v1, const DESreal4vector& v2)
{
double product;
product = v1(0)*v2(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v1(mu)*v2(mu); // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Square of one 4-vector: square(v) = dot(v,v) 
double square(const DESreal4vector& v)
{
double product;
product = v(0)*v(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v(mu)*v(mu); // We use public operator here.
  }
return product;
}
// End of implementation code for DESreal4vector.
//==================================================================================================
// Implementation code for class DEScmplx4vector.

// Default constructor
DEScmplx4vector::DEScmplx4vector()
{
for(int mu = 0; mu <= 3; mu++)
  {
  vector[mu] = 0.0;
  }
}
//--------------------------------------------------------------------------------------------------
// v.in(mu) for input, v.in(mu) = r
complex<double>&  DEScmplx4vector::in(const int& mu)
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// v(mu) for output,  r = v(mu)
const complex<double>&  DEScmplx4vector::operator()(const int& mu) const
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// Sum of two vectors: v3 = v1 + v2
DEScmplx4vector DEScmplx4vector::operator+(const DEScmplx4vector& v2) const
{
DEScmplx4vector sum;
for(int mu = 0; mu <= 3; mu++)
  {
  sum.vector[mu] =  vector[mu] + v2.vector[mu];
  }
return sum;
}
//--------------------------------------------------------------------------------------------------
// Difference of two vectors: v3 = v1 - v2
DEScmplx4vector DEScmplx4vector::operator-(const DEScmplx4vector& v2) const
{
DEScmplx4vector difference;
for(int mu = 0; mu <= 3; mu++)
  {
  difference.vector[mu] =  vector[mu] - v2.vector[mu];
  }
return difference;
}
//--------------------------------------------------------------------------------------------------
// Vector times scalar: v2 = v1*c
DEScmplx4vector  DEScmplx4vector::operator*(const complex<double>& a) const
{
DEScmplx4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.vector[mu] = a * vector[mu];
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// -1 times vector: v2 = - v1
DEScmplx4vector  DEScmplx4vector::operator-() const
{
DEScmplx4vector negative;
for(int mu = 0; mu <= 3; mu++)
  {
  negative.vector[mu] = - vector[mu];
  }
return negative;
}
//--------------------------------------------------------------------------------------------------
// Scalar times vector: v2 = c*v1
DEScmplx4vector operator*(const complex<double>& a, const DEScmplx4vector& v)
{
DEScmplx4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.in(mu) = a * v(mu);  // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Dot product of two 4-vectors: dot(v1,v2)
complex<double> dot(const DEScmplx4vector& v1, const DEScmplx4vector& v2)
{
complex<double> product;
product = v1(0)*v2(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v1(mu)*v2(mu); // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Square of one 4-vector: square(v) = dot(v,v) 
complex<double> square(const DEScmplx4vector& v)
{
complex<double> product;
product = v(0)*v(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v(mu)*v(mu); // We use public operator here.
  }
return product;
}

// End of implementation code for DEScmplx4vector.
//==================================================================================================
// Implementation code for class DESxarray.

// Default constructor
DESxarray::DESxarray()
{
for(int i = 0; i <= maxsize; i++) x[i] = 0.0;
}
//--------------------------------------------------------------------------------------------------
// x.in(i) for input, x.in(i) = ... .
double& DESxarray::in(const int& i)
{
return x[i];
}
//--------------------------------------------------------------------------------------------------
// x(i) for output, ... = x(i).
const double& DESxarray::operator()(const int& i) const
{
return x[i];
}
// End of implementation code for DESxarray.
//==================================================================================================
// Implementation code for class DESzarray.

// Default constructor
DESzarray::DESzarray()
{
for(int i = 0; i <= maxsize; i++) z[i] = 0.0;
}
//--------------------------------------------------------------------------------------------------
// z.in(i) for input, z.in(i) = ... .
complex<double>& DESzarray::in(const int& i)
{
return z[i];
}
//--------------------------------------------------------------------------------------------------
// z(i) for output, ... = z(i).
const complex<double>& DESzarray::operator()(const int& i) const
{
return z[i];
}
// End of implementation code for DESzarray.
//==================================================================================================
// Implementation code for class DES_Smatrix.
// Default constructor
DES_Smatrix::DES_Smatrix()
{
int i,j;
for(i = 0; i < maxsize; i++)
  {
  for(j = 0; j < maxsize; j++) Smatrix[i][j] = 0.0;
  }
}
//--------------------------------------------------------------------------------------------------
// Makes S of type needed.
// Note that we have revised the definition of S_{i,j} so that S_{0,i} = i m_0^2. However
// our S matrix here is real, so we set S[0][i] = 0. Thus we should get m_0^2 from 
// S.getm0sq().
//
void DES_Smatrix::make(DESreal4vector Q[maxsize])
{
double temp;
double tiny = 1.0e-10;
int i,j;
for(i = 0; i <= nverts; i++)
  {
  for(j = 0; j <= nverts; j++) Smatrix[i][j] = 0.0;  // Start with Smatrix set to zero.
  }
//
if(type == PLAIN)
  {
  Smatrix[0][0] = 0.0;
  for(i = 1; i <= nverts; i++)
    {
    Smatrix[0][i] = 0.0;
    Smatrix[i][0] = 0.0;
    Smatrix[i][i] = - twomsq;
    for(j = i+1; j <= nverts; j++)
      {
      temp = square(Q[i] - Q[j]);
      if(abs(temp) < tiny) temp = 0.0;  // correct for roundoff errors
      temp = temp - twomsq;
      Smatrix[i][j] = temp;
      Smatrix[j][i] = temp;
      }
    }
  }
else  // PLAIN is the only type we have for now.
  {
  cout<<"Oops, not programmed for that type of matrix S."<<endl;
  exit(1);
  }
}                              
//--------------------------------------------------------------------------------------------------
// s.in(i,j) for input, s.in(i,j) = x. This is probably not really needed.
double&  DES_Smatrix::in(const int& i, const int& j)
{
return Smatrix[i][j];
}
//--------------------------------------------------------------------------------------------------
// s(i,j) for output,  x = s(i,j).
const double&  DES_Smatrix::operator()(const int& i, const int& j) const
{
return Smatrix[i][j];
}
//--------------------------------------------------------------------------------------------------
// Set rtssq = rts^2.
void DES_Smatrix::setrtssq(double rtssqIN)
{
rtssq = rtssqIN;
}
//--------------------------------------------------------------------------------------------------
// Return rtssq = rts^2.
const double& DES_Smatrix::getrtssq()
{
return rtssq;
}
//--------------------------------------------------------------------------------------------------
// Set rtssq = rts^2.
void DES_Smatrix::setm0sq(double m0sqIN)
{
m0sq = m0sqIN;
}
//--------------------------------------------------------------------------------------------------
// Return m0sq.
const double& DES_Smatrix::getm0sq()
{
return m0sq;
}
//--------------------------------------------------------------------------------------------------
// Set rtssq = rts^2.
void DES_Smatrix::setmass(double massIN)
{
twomsq = 2.0*massIN*massIN;
mass = massIN;
}
//--------------------------------------------------------------------------------------------------
// Return m0sq.
const double& DES_Smatrix::getmass()
{
return mass;
}
//--------------------------------------------------------------------------------------------------
// Set type of term S is used for (PLAIN or UV).
void DES_Smatrix::settype(SubtrTypes typeIN)
{
type = typeIN;
}
//--------------------------------------------------------------------------------------------------
// Return type of term S is used for (PLAIN or UV).
const SubtrTypes& DES_Smatrix::gettype()
{
return type;
}
//--------------------------------------------------------------------------------------------------
void DES_Smatrix::setnverts(int nvertsIN)
{
nverts = nvertsIN;
}
//--------------------------------------------------------------------------------------------------
// Return nverts.
const int& DES_Smatrix::getnverts()
{
return nverts;
}
//
// End of implementation code for DES_Smatrix.
//==================================================================================================
// Implementation code for class DESrandom.
// input : seed (0<=seed<259200)                                            
// output: random numbers (between 2^(-32) and 1-2^(-32) ) 
// Code by D. Sung and DES.

DESrandom::DESrandom(/*in*/int seed) //constructor(initializes seed)   
{
  reSetSeed(seed);
}      
//--------------------------------------------------------------------------------------------------
DESrandom::DESrandom()  //default constructor. (puts seed as 0)
{
  reSetSeed(0);
}      
//--------------------------------------------------------------------------------------------------
bool DESrandom::reSetSeed(/*in*/int seed) 
              //preCondition 0<=seed<259200
              //postCond. : seed is reset, return 0 if seed is invalid.
{
  lbit=31;
  fac=pow(2.0,-lbit);

  int i, ifac1, isk, idi, i1s;
  int im=259200;
  int ia=7141;
  int ic=54773;

  assert(seed >=0);
  assert(seed < im);

  for(int i=0; i<10 ; i++)
    seed=(seed*ia + ic)%im;
  ifac1=((int(pow(2.0,lbit-1))-1)*2+1)/im;
  for(int i=1;i<250;i++)
    {
    seed=(seed*ia+ic)%im;
    ir[i]=ifac1*seed;
    }

  ir[0]=1;
  i=0; 
  isk=250/lbit;
  idi=1;
  i1s=1;

  for(int lb=1;lb<lbit;lb++)
  {
    i=i+isk;
    idi=idi*2;
    i1s=i1s+idi;
    ir[i]=ir[i]|idi;        //bitwise operator. | : or
    ir[i]=ir[i]&i1s;        // & : and
  }

  newran();
  return true;
}
//--------------------------------------------------------------------------------------------------
const double& DESrandom::random()   //return random value.
{
  if(j>=249)
    newran();
  j=j+1;
  return rr[j];
}
//--------------------------------------------------------------------------------------------------   
void DESrandom::newran()
{
  for(int n=0;n<103;n++)
  {
  ir[n]=ir[n+147]^ir[n];
  rr[n]=fac*(double(ir[n])+0.5);
  }
  for(int n=103;n<250;n++)
  {
  ir[n]=ir[n-103]^ir[n];       //^ : bitwise op. exclusive or.
  rr[n]=fac*(double(ir[n])+0.5);
  }

j=-1;
}
// End of implementation code for DESrandom.
//==================================================================================================
// Function to sort an array of real numbers.
void DESsort(double x[], const int& nmax) // Returns sorted array x = {x[0],...,x[nmax-1]}.
{
const int mediumsize = 32;  // (mediumsize = 8 is 20% faster for large nmax). 
int i,j;
double xsmall;
if(nmax < mediumsize)  // action for a small array
  {
  for(i = 0; i < nmax - 1; i++)
    {
    xsmall = x[i];
    for(j = i + 1; j < nmax; j++)
      {
      if(x[j] < xsmall)
        {
        xsmall = x[j];
        x[j] = x[i];
        x[i] = xsmall;
        }  // end if
      }  // end for(j)
    }  // end for(i)
  }
else  // nmax >= mediumsize, so divide x into two parts x1 and x2.
  {
  int n1,n2;
  n1 = nmax/2; // integer divide!
  n2 = nmax - n1;
  double x1[n1];
  double x2[n2];
  for(i = 0; i < n1; i++)
    {
    x1[i] = x[i];
    x2[i] = x[n1+i];
    }
  x2[n2-1] = x[nmax-1];
  DESsort(x1,n1);  // Sort x1.
  DESsort(x2,n2);  // Sort x2.
  // Now merge x1 and x2 into x;\.
  int j1 = 0;
  int j2 = 0;
  for(i = 0; i < nmax; i++)
    {
    if(j1 >= n1)  // We have used up x1, so fill from x2.
      {
      x[i] = x2[j2];
      j2++;
      }
    else if(j2 >= n2)  // We have used up x2, so fill from x1.
      {
      x[i] = x1[j1];
      j1++;
      }
    else if(x1[j1] < x2[j2])  // Select smallest from x1 and x2.
      {
      x[i] = x1[j1];
      j1 ++;
      }
    else
      {
      x[i] = x2[j2];
      j2 ++;
      }  // end if
    }  // end for
  } // end else nmax >= mediumsize so divide x into two parts x1 and x2.
}
//==================================================================================================
// Function n!.
int factorial( const int& n)
{
int result = 1;
for(int i = 2; i <= n; i++ )
  result = result * i;
return result;
}
//==================================================================================================
// Function to give cyclic index in the range 1,...,n.
int cyclic(const int& i, const int& n)
{
int iout;
if(i > 0)
  {
  if(i <= n)
    {
    iout = i;   // 0 < i <= n, the most common case
    }
  else if(i <= 2*n)
    {
    iout = i - n;  // n < i <= 2*n
    }
  else
    {
    iout = (i-1)%n + 1; // 2*n < i , use the general formula
    }
  }
else if(i > -n)
  {
  iout = i + n;  // -n < i <= 0
  }
else
  {
  iout = (i-1)%n + 1; // i <= -n, use the general formula
  }
return iout;
  }
//==================================================================================================
// Function Lambda^2(z).
complex<double> Lambdasq(const int nverts, DES_Smatrix s, DESzarray z)
{
static const complex<double>ii(0.0,1.0);
static int i,j;
static complex<double> Lamsq;
Lamsq = 0;
for(i = 1; i <= nverts; i++)
  {
  for(j = 1; j <= nverts; j++)
    {
    Lamsq = Lamsq + s(i,j)*z(i)*z(j);
    }
  }
Lamsq = 0.5*Lamsq;
Lamsq = Lamsq + ii*z(0)*(1.0 - z(0))*s.getm0sq();
return Lamsq;
}
//==================================================================================================
// Function Lambda^2_abs(x).
// Lambda^2 for real arguments x, calculated with the absolute value of the s(i,j) matrix.
double Lambdasqabs(const int nverts, DES_Smatrix s, DESxarray x)
{
static int i,j;
static double Lamsq;
Lamsq = 0;
for(i = 1; i <= nverts; i++)
  {
  for(j = 1; j <= nverts; j++)
    {
    Lamsq = Lamsq + abs(s(i,j))*x(i)*x(j);
    }
  }
Lamsq = 0.5*Lamsq;
Lamsq = Lamsq + s.getm0sq()*x(0)*(1.0 - x(0));
return Lamsq;
}
//==================================================================================================
// Function Determinant(M,d).
// This function calculates the determinant of any input matrix using LU decomposition.
// Code by W. Gong and D.E. Soper
complex<double> Determinant(const complex<double> bb[maxsize][maxsize], const int& dimension)
{
// Define matrix related variables.
static const complex<double> one(1.0, 0.0);
static complex<double> determinant;
complex<double> aa[dimension][dimension];
int indx[dimension], d;
// Initialize the determinant.
determinant = one;
// Inintialize the matrix to be decomposed with the transferred matrix b.
for(int i=0; i < dimension; i++)
  {
  for(int j=0; j < dimension; j++)
    {
    aa[i][j] = bb[i][j];
    }
  }
// Define parameters used in decomposition.
static int i, imax, j, k, flag=1;
static complex <double> dumc,sum;
static double aamax, dumr;
static double vv[maxsize];
// Initialize the parity parameter.
d = 1;
// Get the implicit scaling information.
for(i = 0; i < dimension; i++)
  {
  aamax=0;
  for(j = 0; j < dimension; j++)
    {
    if(abs(aa[i][j]) > aamax) aamax=abs(aa[i][j]);
    }
    // Set a flag to check if the determinant is zero.
    if(aamax == 0) flag=0;
    // Save the scaling.
    vv[i]=1.0/aamax;
    }
if(flag == 1)
  {
  for(j = 0; j < dimension; j++)
    {
    for(i = 0; i < j; i++)
      {
      sum = aa[i][j];
      for(k = 0; k < i; k++) sum = sum - aa[i][k]*aa[k][j];
      aa[i][j] = sum;
      }
    //Initialize for the search for largest pivot element.
    aamax=0;
    for(i = j; i < dimension; i++)
      {
      sum = aa[i][j];
      for(k = 0; k < j; k++) sum = sum - aa[i][k]*aa[k][j];
        aa[i][j]=sum;
        // Figure of merit for the pivot.
        dumr = vv[i] * abs(sum);
        // Is it better than the best so far?
        if(dumr >= aamax)
          {
          imax=i;
          aamax=dumr;
          }
        }  // End for(i = j; i < dimension; i++)
      // See if we need to interchange rows.
      if(j != imax)
        {
        for(k = 0; k < dimension; k++)
          {
          dumc = aa[imax][k];
          aa[imax][k] = aa[j][k];
          aa[j][k] = dumc;
          }
        // Change the parity of d.
        d = -d;
        // Interchange the scale factor.
        vv[imax] = vv[j];
        }  // End if(j != imax)
      indx[j]=imax;
      if(j != dimension - 1)
        {
        dumc = 1.0/aa[j][j];
        for(i = j+1; i < dimension; i++) aa[i][j] = aa[i][j]*dumc;
        }
    } // End for(j = 0; j < dimension; j++)
  } // END if(flag == 1)
// Calculate the determinant using the decomposed matrix.
if(flag == 0) determinant = 0.0;
else
  {    // Multiply the diagonal elements.
  for(int diagonal = 0; diagonal < dimension; diagonal++)
    {
    determinant = determinant * aa[diagonal][diagonal];
    }
  determinant=(double)d*determinant;
  } // End if(flag == 0) ... else
return determinant;
}
//==================================================================================================
// Function lmod(linfo).
// Transformed loop momentum for numerator function.
DEScmplx4vector lmod(const DES_linfo& linfo)
{
static int i,mu;
static DEScmplx4vector l;          // The loop momentum after all transformations.
static DEScmplx4vector ell;        // The starting vector
static DEScmplx4vector parityell;  // Parity transform of ell.
static const complex<double>ii(0.0,1.0);
static const complex<double>rtmi(0.7071067811865475,-0.7071067811865475);
static complex<double> Lambda;
static complex<double> temp;
//
temp = sqrt(-ii*linfo.Lambdasq*pow(1.0 + ii*linfo.etasum,2))/rtmi;
if(real(temp) < 0.0)
  {
  if (linfo.type == PLAIN)
    {
    cout << "real Lambda*(1.0 + ii*etasum) is negative " << temp << endl;
    exit(1);
    }
  }
if(imag(temp) < 0.0)
  {
  if (linfo.type == PLAIN)
    {
    cout << "imag Lambda*(1.0 + ii*etasum) is negative " << temp << endl;
    exit(1);
    }
  }
Lambda = temp/(1.0 + ii*linfo.etasum);
ell.in(0) = linfo.ell(0);
parityell.in(0) =  linfo.ell(0);
for(mu = 1; mu <= 3; mu++)
  {
  ell.in(mu) = linfo.ell(mu);
  parityell.in(mu) = - linfo.ell(mu);
  }
//
if(linfo.type == PLAIN) 
  {
  l = (1.0 - ii)*ell + (1.0 + ii)*parityell;
  l = 0.5*Lambda*l;
  for(i = 1; i <= linfo.nverts; i++) l = l + linfo.z(i)*linfo.Q[i];
  l = (1.0/(1.0 - linfo.z(0)))*l;
  }
else if(linfo.type == UV)
  {
  l = (1.0 - ii)*ell + (1.0 + ii)*parityell;
  l = 0.5*Lambda*l;
  l = (1.0/(1.0 - linfo.z(0)))*l;
  }
else
 {
  cout<<"Oops, not programmed for that one"<<endl;
  exit(1);
  }
return l;
}
//==================================================================================================
// Dummy numerator function N(ell,Q).
complex<double> numerator(const DEScmplx4vector& lvec, const DEScmplx4vector Q[])
{
static complex<double> thenumerator;
thenumerator = 1.0;
return thenumerator;
}
//==================================================================================================
// Implementation code for class FermionLoop. Authors: W. Gong & D. E. Soper
//
FermionLoop::FermionLoop(int nvertsIN) // Constructor. Initializes gamma.
{
static int i,j;
//
nverts = nvertsIN;
// Initialize polarization vectors to zero. They will be changed later.
for(i = 1; i <= maxsize; i++)
  {
  for(j = 0; j <= 3; j++)
    {
    epsln[i].in(j) = 0;
    }
  }
}
//--------------------------------------------------------------------------------------------------
void FermionLoop::FeedEpsilon(const int nv, const DEScmplx4vector& epsilon)
               // Precondition:
               //    nv is the index for the target vertex.
               //    epsilon is a complex 4-vector.
               // Postcondition:
               //    The polarization vector at the given vertex will be set
               //    equal to the given "epsilon".
{
epsln[nv] = epsilon;
}
//--------------------------------------------------------------------------------------------------
void FermionLoop::CalEpsilonNull(const int nv, const int hIN, const DESreal4vector& pIN)
               // Precondition:
               //    nv is the index for the target vertex.
               //    h could only be either "+1" or "-1", which separately stands
               //    for plus helicity and minus helicity.
               //    l is a 4-vector, q is an array of 4-vectors, where the
               //    external momemtum p[i] = Q[i+1] - Q[i], and loop momemtum
               //    are l-Q[i].
               // Postcondition:
               //    Calculate the polarization vector at the given vertex
               //    with the chosen helicity.
               // This version has null plane gauge.
{
static DESreal4vector p;
static DEScmplx4vector epsilon;
static const complex<double>ii(0.0,1.0);
static complex<double> e1,e2,eps0;
static int h;

p = pIN;
if(p(0) < 0) p = - p; // This gives the physical momentum for the incoming particles.
h = - hIN;  // Gives epsilon^* for outgoing particles and - h convention for incoming particles.
//
if(h == 1)
  {
  e1 =  1.0/sqrt(2.0);   // plus helicity
  e2 =   ii/sqrt(2.0);   // plus helicity
  }
else if(h == -1)
  {
  e1 = 1.0/sqrt(2.0);   // minus helicity
  e2 = - ii/sqrt(2.0);  // minus helicity
  }
else if (h == 0) // return \epsilon^\mu \propto p^\mu.
  {
  epsilon.in(0) = 1.0;
  epsilon.in(1) = p(1)/p(0);
  epsilon.in(2) = p(2)/p(0);
  epsilon.in(3) = p(3)/p(0);
  epsln[nv] = epsilon;
  return;
  }
else
 {
 cout << "Wrong helicity input!" << endl;
 exit(1);
 }
if(p(3) > 0.0) // Use epsilon^+ = 0 vectors.
  {
  eps0 = (p(1)*e1 + p(2)*e2)/(p(0) + p(3));
  epsilon.in(0) = eps0;
  epsilon.in(1) = e1;
  epsilon.in(2) = e2;
  epsilon.in(3) = - eps0;
  }
else                // Use epsilon^+ = 0 vectors.
  {
  e2 = - e2;        // eps_T(h) -> - eps_T(-h).
  eps0 = (p(1)*e1 + p(2)*e2)/(p(0) - p(3));
  epsilon.in(0) = eps0;
  epsilon.in(1) = e1;
  epsilon.in(2) = e2;
  epsilon.in(3) = eps0;
  }
epsln[nv] = epsilon;
}
//--------------------------------------------------------------------------------------------------
void FermionLoop::CalEpsilonCoulomb(const int nv, const int hIN, const DESreal4vector& pIN)
               // Precondition:
               //    nv is the index for the target vertex.
               //    h could only be either "+1" or "-1", which separately stands
               //    for plus helicity and minus helicity.
               //    l is a 4-vector, q is an array of 4-vectors, where the
               //    external momemtum p[i] = Q[i+1] - Q[i], and loop momemtum
               //    are l-Q[i].
               // Postcondition:
               //    Calculate the polarization vector at the given vertex
               //    with the chosen helicity.
               // This version has Coulomb gauge.
{
static DESreal4vector p;
static DEScmplx4vector epsilon;
static double E;
static double r;
static double temp;
static const complex<double>ii(0.0,1.0);
static complex<double> e1,e2;
static int h;

p = pIN;
if(p(0) < 0) p = - p; // This gives the physical momentum for the incoming particles.
h = - hIN;  // This gives epsilon^* for outgoing particles and - h convention for incoming particles.

if(h == 1)
  {
  e1 = 1.0/sqrt(2.0);   // plus helicity
  e2 = ii/sqrt(2.0);   // plus helicity
  }
else if(h == -1)
  {
  e1 = 1.0/sqrt(2.0);   // minus helicity
  e2 = - ii/sqrt(2.0);  // minus helicity
  }
else if (h == 0) // return \epsilon^\mu \propto p^\mu.
  {
  epsilon.in(0) = 1.0;
  epsilon.in(1) = p(1)/p(0);
  epsilon.in(2) = p(2)/p(0);
  epsilon.in(3) = p(3)/p(0);
  epsln[nv] = epsilon;
  return;
  }
else
 {
 cout << "Wrong helicity input!" << endl;
 exit(1);
 }

epsilon.in(0) = 0.0;
epsilon.in(1) = e1;
epsilon.in(2) = e2;
epsilon.in(3) = 0.0;
E = sqrt( p(1)*p(1) + p(2)*p(2) + p(3)*p(3) );
r = sqrt( p(1)*p(1) + p(2)*p(2) );
if(r/E > 1.0e-10)
  {
  epsilon.in(1) = p(2)*(e1*p(2)-e2*p(1))/(r*r) + p(1)*p(3)*(e1*p(1)+e2*p(2))/(E*r*r);
  epsilon.in(2) = -p(1)*(e1*p(2)-e2*p(1))/(r*r) + p(2)*p(3)*(e1*p(1)+e2*p(2))/(E*r*r);
  temp = 1.0 - pow((p(3)/E), 2);
  epsilon.in(3) = -sqrt(temp)*(e1*p(1)+e2*p(2))/r;
  }
else if(p(3) < 0)        // Momentum along the negative z-axis.
  {
  epsilon.in(1) = - e1;  // eps_T(h) -> - eps_T(-h).
  epsilon.in(2) =   e2;
  }
epsln[nv] = epsilon;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorV(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static complex<double> ii(0.0,1.0);
static complex<double> numerator;
static int nv,i,j,mu;
static complex<double> eslash[5][5];
static complex<double> kslash[5][5];
static complex<double> product[5][5];
static complex<double> product1[5][5];
static DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// 
// Start with the unit matrix.
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 1; i <= 4; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  kslash[1][3] = k(0) - k(3);
  kslash[1][4] = -k(1) + ii*k(2);
  kslash[2][3] = -k(1) - ii*k(2);
  kslash[2][4] = k(0) + k(3);
  //
  kslash[3][1] = k(0) + k(3);
  kslash[3][2] = k(1) - ii*k(2);
  kslash[4][1] = k(1) + ii*k(2);
  kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash.
  product1[1][3] = kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  product1[1][4] = kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  product1[2][3] = kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  product1[2][4] = kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1];
  product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2];
  product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1];
  product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2];
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  eslash[1][3] = eps(0) - eps(3);
  eslash[1][4] = -eps(1) + ii*eps(2);
  eslash[2][3] = -eps(1) - ii*eps(2);
  eslash[2][4] = eps(0) + eps(3);
  //
  eslash[3][1] = eps(0) + eps(3);
  eslash[3][2] = eps(1) - ii*eps(2);
  eslash[4][1] = eps(1) + ii*eps(2);
  eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  //
  product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace.
numerator = 0.0;
for(i = 1; i <= 4; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorM(const DEScmplx4vector& l, const DEScmplx4vector Q[], const double m)
{
static complex<double> ii(0.0,1.0);
static complex<double> numerator;
static int nv,i,j,mu;
static complex<double> eslash[5][5];
static complex<double> kslash[5][5];
static complex<double> product[5][5];
static complex<double> product1[5][5];
static DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators, with mass:
// Trace{ epsilonslash[nverts]*(kslash[nverts] + m)*...*epsilonslash[1]*(kslash[1] + m) }
// 
// Start with the unit matrix.
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 1; i <= 4; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  kslash[1][3] = k(0) - k(3);
  kslash[1][4] = -k(1) + ii*k(2);
  kslash[2][3] = -k(1) - ii*k(2);
  kslash[2][4] = k(0) + k(3);
  //
  kslash[3][1] = k(0) + k(3);
  kslash[3][2] = k(1) - ii*k(2);
  kslash[4][1] = k(1) + ii*k(2);
  kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash + m.
  product1[1][1] = m*product[1][1] + kslash[1][3]*product[3][1] + kslash[1][4]*product[4][1];
  product1[1][2] = m*product[1][2] + kslash[1][3]*product[3][2] + kslash[1][4]*product[4][2];
  product1[1][3] = m*product[1][3] + kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  product1[1][4] = m*product[1][4] + kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  //
  product1[2][1] = m*product[2][1] + kslash[2][3]*product[3][1] + kslash[2][4]*product[4][1];
  product1[2][2] = m*product[2][2] + kslash[2][3]*product[3][2] + kslash[2][4]*product[4][2];
  product1[2][3] = m*product[2][3] + kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  product1[2][4] = m*product[2][4] + kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1] + m*product[3][1];
  product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2] + m*product[3][2];
  product1[3][3] = kslash[3][1]*product[1][3] + kslash[3][2]*product[2][3] + m*product[3][3];
  product1[3][4] = kslash[3][1]*product[1][4] + kslash[3][2]*product[2][4] + m*product[3][4];
//
  product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1] + m*product[4][1];
  product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2] + m*product[4][2];
  product1[4][3] = kslash[4][1]*product[1][3] + kslash[4][2]*product[2][3] + m*product[4][3];
  product1[4][4] = kslash[4][1]*product[1][4] + kslash[4][2]*product[2][4] + m*product[4][4];
  //
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  eslash[1][3] = eps(0) - eps(3);
  eslash[1][4] = -eps(1) + ii*eps(2);
  eslash[2][3] = -eps(1) - ii*eps(2);
  eslash[2][4] = eps(0) + eps(3);
  //
  eslash[3][1] = eps(0) + eps(3);
  eslash[3][2] = eps(1) - ii*eps(2);
  eslash[4][1] = eps(1) + ii*eps(2);
  eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  product[1][3] = eslash[1][3]*product1[3][3] + eslash[1][4]*product1[4][3];
  product[1][4] = eslash[1][3]*product1[3][4] + eslash[1][4]*product1[4][4];
  //
  product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  product[2][3] = eslash[2][3]*product1[3][3] + eslash[2][4]*product1[4][3];
  product[2][4] = eslash[2][3]*product1[3][4] + eslash[2][4]*product1[4][4];
  //
  product[3][1] = eslash[3][1]*product1[1][1] + eslash[3][2]*product1[2][1];
  product[3][2] = eslash[3][1]*product1[1][2] + eslash[3][2]*product1[2][2];
  product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  //
  product[4][1] = eslash[4][1]*product1[1][1] + eslash[4][2]*product1[2][1];
  product[4][2] = eslash[4][1]*product1[1][2] + eslash[4][2]*product1[2][2];
  product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace.
numerator = 0.0;
for(i = 1; i <= 4; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorUV(const DEScmplx4vector& l)
{
// This one if for the "extra" part of the UV subtraction.
// This is only for nverts = 4 with vector vertices.
static int i;
static complex<double> numerator;
numerator = 32.0;
for(i = 1; i <= 4; i++) numerator = numerator * dot(epsln[i],l);
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::AltNumeratorV(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static DEScmplx4vector v[2*maxsize];
static int i;
static complex<double> numerator;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// In this case, we expand the trace in dot products of 4-vectors.
//
for(i = 1; i<= nverts; i++)
  {
  v[2*(nverts - i + 1)]     = l - Q[i];
  v[2*(nverts - i + 1) - 1] = epsln[i];
  }
numerator = DiracTrace(v,2*nverts); // Expands the trace in dot products of 4 vectors.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorL(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static complex<double> ii(0.0,1.0);
static complex<double> numerator;
static int nv,i,j,mu;
static complex<double> eslash[5][5];
static complex<double> kslash[5][5];
static complex<double> product[5][5];
static complex<double> product1[5][5];
static DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// 
// Start with the matrix diag(0,0,1,1).
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 3; i <= 4; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  kslash[1][3] = k(0) - k(3);
  kslash[1][4] = -k(1) + ii*k(2);
  kslash[2][3] = -k(1) - ii*k(2);
  kslash[2][4] = k(0) + k(3);
  // We don't need the following part of kslash.
  //kslash[3][1] = k(0) + k(3);
  //kslash[3][2] = k(1) - ii*k(2);
  //kslash[4][1] = k(1) + ii*k(2);
  //kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash.
  product1[1][3] = kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  product1[1][4] = kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  product1[2][3] = kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  product1[2][4] = kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  //product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1];
  //product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2];
  //product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1];
  //product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2];
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  // The following part of epsilon slash is removed by the left-handed projection.
  //eslash[1][3] = eps(0) - eps(3);
  //eslash[1][4] = -eps(1) + ii*eps(2);
  //eslash[2][3] = -eps(1) - ii*eps(2);
  //eslash[2][4] = eps(0) + eps(3);
  //
  eslash[3][1] = eps(0) + eps(3);
  eslash[3][2] = eps(1) - ii*eps(2);
  eslash[4][1] = eps(1) + ii*eps(2);
  eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  //product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  //product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  //product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  //product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  //
  product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace of the {3,4} diagonal elements.
numerator = 0.0;
for(i = 3; i <= 4; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorR(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static complex<double> ii(0.0,1.0);
static complex<double> numerator;
static int nv,i,j,mu;
static complex<double> eslash[5][5];
static complex<double> kslash[5][5];
static complex<double> product[5][5];
static complex<double> product1[5][5];
static DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// 
// Start with the matrix diag(1,1,0,0).
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 1; i <= 2; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  //  We don't need the following part of kslash.
  //kslash[1][3] = k(0) - k(3);
  //kslash[1][4] = -k(1) + ii*k(2);
  //kslash[2][3] = -k(1) - ii*k(2);
  //kslash[2][4] = k(0) + k(3);
  //
  kslash[3][1] = k(0) + k(3);
  kslash[3][2] = k(1) - ii*k(2);
  kslash[4][1] = k(1) + ii*k(2);
  kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash.
  //product1[1][3] = kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  //product1[1][4] = kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  //product1[2][3] = kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  //product1[2][4] = kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1];
  product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2];
  product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1];
  product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2];
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  //
  eslash[1][3] = eps(0) - eps(3);
  eslash[1][4] = -eps(1) + ii*eps(2);
  eslash[2][3] = -eps(1) - ii*eps(2);
  eslash[2][4] = eps(0) + eps(3);
  // The following part of epsilon slash is removed by the right-handed projection.
  //eslash[3][1] = eps(0) + eps(3);
  //eslash[3][2] = eps(1) - ii*eps(2);
  //eslash[4][1] = eps(1) + ii*eps(2);
  //eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  //
  //product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  //product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  //product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  //product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace of the {1,2} diagonal elements.
numerator = 0.0;
for(i = 1; i <= 2; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
// End of implementation code for class FermionLoop.
//==================================================================================================
// Function to generate the trace of v1slash*v2slash*...*vnslash.
complex<double> DiracTrace(DEScmplx4vector v[2*maxsize], int n)
{
static complex<double> ii(0.0,1.0);
complex<double> thetrace;
int i,j;
DEScmplx4vector u[10];
complex<double> sign;
//
if(n == 1) thetrace = 0.0;
else if(n == 2) thetrace = 4.0*dot(v[1],v[2]);
else
{
for(i = 1; i <= n - 2; i++) u[i] = v[i+2];
sign = 1.0;
thetrace = 0.0;
for(i = 2; i <= n; i++)
  {
  thetrace = thetrace + sign*dot(v[1],v[i])*DiracTrace(u, n-2);
  for(j = 1; j <= n - 3; j++) u[j] = u[j + 1];
  u[n-2] = v[i];
  sign = - sign;
  }
}
return thetrace;
}
//==================================================================================================
// Function getperms(n,permlist).
// Generates a list of permutations of the integers 1,2,...,n.
// The ith permutation is pi[j] = PermList[i][j].
// For j > n, getperms gives pi[j] = j.
void getperms(int n, int PermList[maxperms][maxsize])
{
static int pi[maxsize];
static int count;
static bool up;
static int i;
if(n > 7)
  {
  cout << "Actually, getperms is not set up to permute more than seven objects." << endl;
  exit(1);
  }
for(i = 1; i <= maxsize; i++) pi[i] = i;
pi[0] = 0;
up = false;
count = 0;
perm(n,pi,count,up,PermList);
}
//--------------
// This calls the recursively defined function perm(...).
void perm(int n, int pi[], int& count, bool up, int PermList[maxperms][maxsize])
{
int i,j,k,temp;
// If called from "higher n" version of perm and here n > 2,
// call perm() to construct permutations of n-1 objects.
if(!up && n > 2) perm(n-1,pi,count,up,PermList);
// Now that we have constructed permuations of (n - 1) objects, we want to
// rotate the first n elements and then construct permutations of new list of
// (n - 1) objects. We will do this for all (n-1) rotated version of our
// n objects (in addition to the unrotated version).
for (i = 1; i < n; i++)
  {
  if (n == 2)
    {
    // Increment count and record this permutation.
    count++;                                
    PermList[count][0] = count;
    for(k = 1; k < maxsize; k++) PermList[count][k] = pi[k]; 
    // Rotate first n elements, here for n = 2.
    temp = pi[1];                            
    pi[1] = pi[2];
    pi[2] = temp;
    count++;
    // Increment count and record this permutation.
    PermList[count][0] = count;
    for(k = 1; k < maxsize; k++) PermList[count][k] = pi[k]; 
    }
  if(n > 2)
    {
    // Rotate first n elements.
    temp = pi[1];
    for (j = 1; j < n; j++) pi[j] = pi[j+1];
    pi[n] = temp;
    // Call perm() to construct permutations of n-1 objects.
    up = false;
    perm(n-1,pi,count,up,PermList);
    }
  }
// We have rotated first n elements (n - 1) times.
// Rotate first n elements once more.
temp = pi[1];
for (j = 1; j < n; j++) pi[j] = pi[j+1];
pi[n] = temp;
up = true;                             // Now we are back at original state.
return;                                // We return "up" to perm with one greater n.
}
// End of function perm(n, pi[], count, up, PermList)
//==================================================================================================
// Function fourphoton(helicities, pvectors).
// Produces the four-photon scattering amplitude.
// Result is from
//   T.~Binoth, E.~W.~N.~Glover, P.~Marquard and J.~J.~van der Bij,
//   %``Two-loop corrections to light-by-light scattering in supersymmetric  QED,''
//   JHEP {\bf 0205}, 060 (2002)
//   [arXiv:hep-ph/0202266].
//
double fourphoton(int helicity[maxsize],DESreal4vector p[maxsize])
{
double s,t,u;
double X,Y;
double calM;
const double pi = 3.14159265358979;
//-------
s = 2.0*dot(p[1],p[2]);
t = 2.0*dot(p[1],p[3]);
u = 2.0*dot(p[1],p[4]);
X = log(-t/s);
Y = log(-u/s);
if(helicity[1] == 1 && helicity[2] == 1 && helicity[3] == 1 && helicity[4] == 1)
  {
  calM = -8.0;
  }
else if(helicity[1] == 1 && helicity[2] == 1 && helicity[3] == 1 && helicity[4] == -1)
  {
  calM = 8.0;
  }
else if(helicity[1] == 1 && helicity[2] == 1 && helicity[3] == -1 && helicity[4] == -1)
  {
  calM =      - 4.0 - 4.0*((X-Y)*(X-Y) + pi*pi - 2*(X-Y))*t*t/s/s;
  calM = calM - 4.0 - 4.0*((Y-X)*(Y-X) + pi*pi - 2*(Y-X))*u*u/s/s;
  }
else
  {
  cout << "Actually, I am not prepared for that helicity combination." << endl;
  calM = 0.0;
  }
return calM;
}
//==================================================================================================