minorex.cpp

Go to the documentation of this file.

/*
 * minorex.cpp - extra integrals for asymptotic expansion
 *
 * this file is part of PJFry library
 * Copyright 2011 Valery Yundin
 */
 
#include "cache.h"
#include "minor.h"
#include <stdio.h>
 
/* ========================================================
 *               IR triangles
 * ========================================================
 */
 
/* --------------------------------------------------------
 *   I2mDstu bubble in D-2 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I2mDstu( int ep, int s, int t, int u, int m, int n ) {
  ncomplex     sum1    = 0;
  const double dstustu = M3( s, t, u, s, t, u );
  const double msq1    = kinem.mass( m );
  const double msq2    = kinem.mass( n );
  if ( ep == 0 )
  {
    const double s_cutoff = seps1 * pmaxM2[im2( m, n ) - 5];
 
    if ( fabs( dstustu ) <= s_cutoff )
    {
      const double mm12 = msq1 - msq2;
      if ( fabs( mm12 ) < meps )
      {
        if ( msq1 > meps ) { sum1 = 1 / msq1; }
        else { sum1 = 0; }
      }
      else
      {
        sum1 = ( -( msq1 > meps ? ICache::getI1( ep, Kinem1( msq1 ) ) / msq1 : 1. ) +
                 ( msq2 > meps ? ICache::getI1( ep, Kinem1( msq2 ) ) / msq2 : 1. ) ) /
               ( mm12 );
      }
    }
    else
    {
      ncomplex sumX = I2stu( ep, s, t, u ) - 2. * I2stu( ep + 1, s, t, u );
      ncomplex sumY = msq2 > meps ? 1. - ICache::getI1( ep, Kinem1( msq2 ) ) / msq2 : 0;
      ncomplex sumZ = msq1 > meps ? 1. - ICache::getI1( ep, Kinem1( msq1 ) ) / msq1 : 0;
 
      const double ds0tu =
          ( Cay[nss( m, m )] * Cay[nss( n, n )] - Cay[ns( m, n )] * Cay[ns( m, n )] );
      sum1 += sumX * dstustu;
 
      const double dsvtuY =
          -( Cay[nss( n, n )] - Cay[ns( m, n )] ); /* minus sign of minor v=m */
      sum1 += sumY * dsvtuY;
 
      const double dsvtuZ =
          +( Cay[ns( m, n )] - Cay[nss( m, m )] ); /* plus sign of minor v=n */
      sum1 += sumZ * dsvtuZ;
 
      sum1 /= ds0tu;
    }
  }
  else if ( ep == 1 )
  {
    if ( fabs( msq1 ) < meps )
    {
      if ( fabs( msq2 ) < meps )
      {
        if ( fabs( dstustu ) > meps )
        {
          sum1 = 2. / ( -0.5 * dstustu ); // I(s,0,0)
        }
        else
        {
          sum1 = 0; // I(0,0,0)
        }
      }
      else
      {
        sum1 = 1. / ( ( -0.5 * dstustu ) - msq2 ); // I(s,0,m2)
      }
    }
    else if ( fabs( msq2 ) < meps )
    {
      sum1 = 1. / ( ( -0.5 * dstustu ) - msq1 ); // I(s,m1,0)
    }
  }
  return sum1;
}
 
/* --------------------------------------------------------
 *   I2stui bubble in D dim with a dot
 * --------------------------------------------------------
 */
ncomplex Minor5::I2stui( int ep, int s, int t, int u, int i, int ip ) {
  ncomplex     sum1    = 0;
  const double dstustu = M3( s, t, u, s, t, u );
  const double msq1    = kinem.mass( i );
  const double msq2    = kinem.mass( ip );
  if ( ep == 0 )
  {
    const double s_cutoff = seps1 * pmaxM2[im2( i, ip ) - 5];
 
    if ( fabs( dstustu ) <= s_cutoff )
    {
      const double mm12 = msq1 - msq2;
      if ( fabs( mm12 ) < meps )
      {
        if ( msq1 > meps ) { sum1 = -1 / ( 2 * msq1 ); }
        else { sum1 = 0; }
      }
      else
      {
        if ( msq1 > meps )
        {
          sum1 = -1 / ( msq1 - msq2 ) - ( +msq2 * ICache::getI1( ep, Kinem1( msq1 ) ) -
                                          msq1 * ICache::getI1( ep, Kinem1( msq2 ) ) ) /
                                            ( msq1 * mm12 * mm12 );
        }
        else { sum1 = ICache::getI1( ep, Kinem1( msq2 ) ) / ( msq2 * msq2 ); }
      }
    }
    else
    {
      sum1 += +( Cay[nss( ip, ip )] - Cay[ns( i, ip )] ) * I2mDstu( ep, s, t, u, i, ip );
      sum1 += msq1 > meps ? 1. - ICache::getI1( ep, Kinem1( msq1 ) ) / msq1 : 0;
      sum1 -= msq2 > meps ? 1. - ICache::getI1( ep, Kinem1( msq2 ) ) / msq2 : 0;
      sum1 /= dstustu;
    }
  }
  else if ( ep == 1 )
  {
    if ( fabs( msq1 ) < meps )
    {
      if ( fabs( dstustu ) > meps )
      {
        if ( fabs( msq2 ) < meps )
        {
          sum1 = -1. / ( -0.5 * dstustu ); // I(s,0,0)
        }
        else
        {
          sum1 = -1. / ( ( -0.5 * dstustu ) - msq2 ); // I(s,0,m2)
        }
      }
      else
      {
        sum1 = 1. / msq2; // I(0,0,m2)
      }
    }
  }
  return sum1;
}
 
/* ========================================================
 *               SMALL Gram4 expansion
 * ========================================================
 */
 
/* --------------------------------------------------------
 *   I2D3stu bubble in D+6 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I2D3stu( int ep, int s, int t, int u ) {
  assert( t != u && u != s && s != t ); // if (t==u || u==s || s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I2D3stu + ep] )
  {
    I2D3stuEval( 0, ep, 1, 2, 3, 4, 5, kinem.p5() );
    I2D3stuEval( 1, ep, 1, 2, 4, 3, 5, kinem.s45() );
    I2D3stuEval( 2, ep, 1, 2, 5, 3, 4, kinem.p4() );
 
    I2D3stuEval( 3, ep, 1, 3, 4, 2, 5, kinem.s12() );
    I2D3stuEval( 4, ep, 1, 3, 5, 2, 4, kinem.s34() );
 
    I2D3stuEval( 5, ep, 1, 4, 5, 2, 3, kinem.p3() );
 
    if ( smax == 5 )
    {
      I2D3stuEval( 6, ep, 2, 3, 4, 1, 5, kinem.p1() );
      I2D3stuEval( 7, ep, 2, 3, 5, 1, 4, kinem.s15() );
 
      I2D3stuEval( 8, ep, 2, 4, 5, 1, 3, kinem.s23() );
 
      I2D3stuEval( 9, ep, 3, 4, 5, 1, 2, kinem.p2() );
    }
 
    fEval[E_I2D3stu + ep] = true;
  }
  int idx = im3( s, t, u ) - 10;
  return pI2D3stu[ep][idx];
}
 
void Minor5::I2D3stuEval( int idx, int ep, int s, int t, int u, int m, int n, double qsq ) {
  ncomplex sum1 = 0;
  if ( ep == 0 )
  {
    const double dstustu  = -2 * qsq; /*M3(s,t,u,s,t,u);*/
    const double msq1     = kinem.mass( m );
    const double msq2     = kinem.mass( n );
    const double s_cutoff = seps1 * pmaxM2[im2( m, n ) - 5];
 
    if ( fabs( dstustu ) <= s_cutoff )
    {
      const double mm12 = msq1 - msq2;
      if ( fabs( mm12 ) < meps )
      { sum1 = -( 5 * msq1 + 6. * ICache::getI1( ep, Kinem1( msq1 ) ) ) * msq1 * msq1 / 36.; }
      else
      {
        sum1 = -13 * ( ( msq1 + msq2 ) * ( msq1 * msq1 + msq2 * msq2 ) ) / 288 +
               ( -msq1 * msq1 * msq1 * ICache::getI1( ep, Kinem1( msq1 ) ) +
                 msq2 * msq2 * msq2 * ICache::getI1( ep, Kinem1( msq2 ) ) ) /
                   ( 24 * mm12 );
      }
    }
    else
    {
      ncomplex sumX = 7. * I2D2stu( ep, s, t, u ) + 2. * I2D2stu( ep + 1, s, t, u );
      ncomplex sumY =
          ( 42. * ICache::getI1( ep, Kinem1( msq2 ) ) + 47 * msq2 ) * msq2 * msq2 / 36.;
      ncomplex sumZ =
          ( 42. * ICache::getI1( ep, Kinem1( msq1 ) ) + 47 * msq1 ) * msq1 * msq1 / 36.;
 
      const double ds0tu =
          ( Cay[nss( m, m )] * Cay[nss( n, n )] - Cay[nss( m, n )] * Cay[nss( m, n )] );
      sum1 += sumX * ds0tu;
 
      const double dsvtuY =
          -( Cay[nss( n, n )] - Cay[nss( m, n )] ); /* minus sign of minor v=m */
      sum1 -= sumY * dsvtuY;
 
      const double dsvtuZ =
          +( Cay[nss( m, n )] - Cay[nss( m, m )] ); /* plus sign of minor v=n */
      sum1 -= sumZ * dsvtuZ;
 
      sum1 /= 49 * dstustu;
    }
  }
  else
  {
    assert( ep == 1 );
    const double y11 = Cay[nss( m, m )];
    const double y12 = Cay[nss( m, n )];
    const double y22 = Cay[nss( n, n )];
    sum1 =
        -( +y11 * y11 * ( y22 + 5 * ( y11 + y12 ) ) + y22 * y22 * ( y11 + 5 * ( y22 + y12 ) ) +
           2 * y12 * y12 * ( y12 + 2 * ( y11 + y22 ) ) + 3 * y11 * y22 * y12 ) /
        1680;
  }
  pI2D3stu[ep][idx] = sum1;
}
 
/* --------------------------------------------------------
 *   I3D4st triangle in D+8 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I3D4st( int ep, int s, int t ) {
  assert( s != t ); // if (s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I3D4st + ep] ) { I3D4stEval( ep ); }
  int idx = im2( s, t ) - 5;
  return pI3D4st[ep][idx];
}
 
void Minor5::I3D4stEval( int ep ) {
  for ( int s = 1; s <= smax; s++ )
  {
    for ( int t = s + 1; t <= 5; t++ )
    {
      int idx = im2( s, t ) - 5;
 
      const double dstst  = M2( s, t, s, t );
      ncomplex     ivalue = 0;
 
      if ( pmaxS3[idx] <= deps )
      {
        printf( "I3D4 - M2({%d,%d}) <= eps\n", s, t );
        ivalue = std::numeric_limits<double>::quiet_NaN();
      }
      else
      {
        double cf = 1. / 8.;
        for ( int ei = ep; ei <= 1; ei++ )
        {
          ncomplex sum1 = M3( 0, s, t, 0, s, t ) * I3D3st( ei, s, t );
          for ( int u = 1; u <= 5; u++ )
          {
            if ( t == u || s == u ) continue;
            sum1 -= M3( u, s, t, 0, s, t ) * I2D3stu( ei, s, t, u );
          }
          ivalue += cf * sum1;
          cf *= 1. / 4.;
        }
        ivalue /= dstst;
      }
      pI3D4st[ep][idx] = ivalue;
    }
  }
  fEval[E_I3D4st + ep] = true;
}
 
/* --------------------------------------------------------
 *   I2D4stu bubble in D+8 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I2D4stu( int ep, int s, int t, int u ) {
  assert( t != u && u != s && s != t ); // if (t==u || u==s || s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I2D4stu + ep] )
  {
    I2D4stuEval( 0, ep, 1, 2, 3, 4, 5, kinem.p5() );
    I2D4stuEval( 1, ep, 1, 2, 4, 3, 5, kinem.s45() );
    I2D4stuEval( 2, ep, 1, 2, 5, 3, 4, kinem.p4() );
 
    I2D4stuEval( 3, ep, 1, 3, 4, 2, 5, kinem.s12() );
    I2D4stuEval( 4, ep, 1, 3, 5, 2, 4, kinem.s34() );
 
    I2D4stuEval( 5, ep, 1, 4, 5, 2, 3, kinem.p3() );
 
    if ( smax == 5 )
    {
      I2D4stuEval( 6, ep, 2, 3, 4, 1, 5, kinem.p1() );
      I2D4stuEval( 7, ep, 2, 3, 5, 1, 4, kinem.s15() );
 
      I2D4stuEval( 8, ep, 2, 4, 5, 1, 3, kinem.s23() );
 
      I2D4stuEval( 9, ep, 3, 4, 5, 1, 2, kinem.p2() );
    }
 
    fEval[E_I2D4stu + ep] = true;
  }
  int idx = im3( s, t, u ) - 10;
  return pI2D4stu[ep][idx];
}
 
void Minor5::I2D4stuEval( int idx, int ep, int s, int t, int u, int m, int n, double qsq ) {
  ncomplex sum1 = 0;
  if ( ep == 0 )
  {
    const double dstustu  = -2 * qsq; /*M3(s,t,u,s,t,u);*/
    const double msq1     = kinem.mass( m );
    const double msq2     = kinem.mass( n );
    const double s_cutoff = seps1 * pmaxM2[im2( m, n ) - 5];
 
    if ( fabs( dstustu ) <= s_cutoff )
    {
      const double mm12 = msq1 - msq2;
      if ( fabs( mm12 ) < meps )
      {
        sum1 = ( 13 * msq1 + 12. * ICache::getI1( ep, Kinem1( msq1 ) ) ) * msq1 * msq1 * msq1 /
               288.;
      }
      else
      {
        const double msq1p4 = ( msq1 * msq1 ) * ( msq1 * msq1 );
        const double msq2p4 = ( msq2 * msq2 ) * ( msq2 * msq2 );
        sum1                = ( 77 * ( msq1 * msq1p4 - msq2 * msq2p4 ) / 7200 +
                 ( +msq1p4 * ICache::getI1( ep, Kinem1( msq1 ) ) -
                   msq2p4 * ICache::getI1( ep, Kinem1( msq2 ) ) ) /
                     120. ) /
               mm12;
      }
    }
    else
    {
      ncomplex sumX = 9. * I2D3stu( ep, s, t, u ) + 2. * I2D3stu( ep + 1, s, t, u );
      ncomplex sumY = -( 36. * ICache::getI1( ep, Kinem1( msq2 ) ) + 47 * msq2 ) * msq2 *
                      msq2 * msq2 / 96.;
      ncomplex sumZ = -( 36. * ICache::getI1( ep, Kinem1( msq1 ) ) + 47 * msq1 ) * msq1 *
                      msq1 * msq1 / 96.;
 
      const double ds0tu =
          ( Cay[nss( m, m )] * Cay[nss( n, n )] - Cay[nss( m, n )] * Cay[nss( m, n )] );
      sum1 += sumX * ds0tu;
 
      const double dsvtuY =
          -( Cay[nss( n, n )] - Cay[nss( m, n )] ); /* minus sign of minor v=m */
      sum1 -= sumY * dsvtuY;
 
      const double dsvtuZ =
          +( Cay[nss( m, n )] - Cay[nss( m, m )] ); /* plus sign of minor v=n */
      sum1 -= sumZ * dsvtuZ;
 
      sum1 /= 81 * dstustu;
    }
    /* printf("I2D4stu(%e,%e,%e)^%d =
     * %e,%e\n",-0.5*dstustu,msq1,msq2,ep,sum1.real(),sum1.imag());
     */
  }
  else
  {
    assert( ep == 1 );
    const double y11 = Cay[nss( m, m )];
    const double y12 = Cay[nss( m, n )];
    const double y22 = Cay[nss( n, n )];
    sum1             = ( +y11 * y11 *
                 ( y11 * ( 35 * ( y11 + y12 ) + 5 * y22 ) + 15 * y12 * ( 2 * y12 + y22 ) ) +
             y22 * y22 *
                 ( y22 * ( 35 * ( y22 + y12 ) + 5 * y11 ) + 15 * y12 * ( 2 * y12 + y11 ) ) +
             y12 * y12 * ( y12 * ( 8 * y12 + 20 * y11 + 20 * y22 ) + 24 * y11 * y22 ) +
             3 * y11 * y11 * y22 * y22 ) /
           120960;
  }
  pI2D4stu[ep][idx] = sum1;
}
 
/* --------------------------------------------------------
 *   I3D5st triangle in D+10 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I3D5st( int ep, int s, int t ) {
  assert( s != t ); // if (s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I3D5st + ep] ) { I3D5stEval( ep ); }
  int idx = im2( s, t ) - 5;
  return pI3D5st[ep][idx];
}
 
void Minor5::I3D5stEval( int ep ) {
  for ( int s = 1; s <= smax; s++ )
  {
    for ( int t = s + 1; t <= 5; t++ )
    {
      int idx = im2( s, t ) - 5;
 
      const double dstst  = M2( s, t, s, t );
      ncomplex     ivalue = 0;
 
      if ( pmaxS3[idx] <= deps )
      {
        printf( "I3D5 - M2({%d,%d}) <= eps\n", s, t );
        ivalue = std::numeric_limits<double>::quiet_NaN();
      }
      else
      {
        double cf = 1. / 10.;
        for ( int ei = ep; ei <= 1; ei++ )
        {
          ncomplex sum1 = M3( 0, s, t, 0, s, t ) * I3D4st( ei, s, t );
          for ( int u = 1; u <= 5; u++ )
          {
            if ( t == u || s == u ) continue;
            sum1 -= M3( u, s, t, 0, s, t ) * I2D4stu( ei, s, t, u );
          }
          ivalue += cf * sum1;
          cf *= 1. / 5.;
        }
        ivalue /= dstst;
      }
      pI3D5st[ep][idx] = ivalue;
    }
  }
  fEval[E_I3D5st + ep] = true;
}
 
/* --------------------------------------------------------
 *   I2D5stu bubble in D+10 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I2D5stu( int ep, int s, int t, int u ) {
  assert( t != u && u != s && s != t ); // if (t==u || u==s || s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I2D5stu + ep] )
  {
    I2D5stuEval( 0, ep, 1, 2, 3, 4, 5, kinem.p5() );
    I2D5stuEval( 1, ep, 1, 2, 4, 3, 5, kinem.s45() );
    I2D5stuEval( 2, ep, 1, 2, 5, 3, 4, kinem.p4() );
 
    I2D5stuEval( 3, ep, 1, 3, 4, 2, 5, kinem.s12() );
    I2D5stuEval( 4, ep, 1, 3, 5, 2, 4, kinem.s34() );
 
    I2D5stuEval( 5, ep, 1, 4, 5, 2, 3, kinem.p3() );
 
    if ( smax == 5 )
    {
      I2D5stuEval( 6, ep, 2, 3, 4, 1, 5, kinem.p1() );
      I2D5stuEval( 7, ep, 2, 3, 5, 1, 4, kinem.s15() );
 
      I2D5stuEval( 8, ep, 2, 4, 5, 1, 3, kinem.s23() );
 
      I2D5stuEval( 9, ep, 3, 4, 5, 1, 2, kinem.p2() );
    }
 
    fEval[E_I2D5stu + ep] = true;
  }
  int idx = im3( s, t, u ) - 10;
  return pI2D5stu[ep][idx];
}
 
void Minor5::I2D5stuEval( int idx, int ep, int s, int t, int u, int m, int n, double qsq ) {
  ncomplex sum1 = 0;
  if ( ep == 0 )
  {
    const double dstustu  = -2 * qsq; /*M3(s,t,u,s,t,u);*/
    const double msq1     = kinem.mass( m );
    const double msq2     = kinem.mass( n );
    const double s_cutoff = seps1 * pmaxM2[im2( m, n ) - 5];
 
    if ( fabs( dstustu ) <= s_cutoff )
    {
      const double mm12 = msq1 - msq2;
      if ( fabs( mm12 ) < meps )
      {
        sum1 = -( 77 * msq1 + 60. * ICache::getI1( ep, Kinem1( msq1 ) ) ) * ( msq1 * msq1 ) *
               ( msq1 * msq1 ) / 7200.;
      }
      else
      {
        const double msq1p5 = ( msq1 * msq1 ) * ( msq1 * msq1 ) * msq1;
        const double msq2p5 = ( msq2 * msq2 ) * ( msq2 * msq2 ) * msq2;
        sum1                = -( 29 * ( msq1 * msq1p5 - msq2 * msq2p5 ) / 14400 +
                  ( +msq1p5 * ICache::getI1( ep, Kinem1( msq1 ) ) -
                    msq2p5 * ICache::getI1( ep, Kinem1( msq2 ) ) ) /
                      720. ) /
               mm12;
      }
    }
    else
    {
      ncomplex sumX = 11. * I2D4stu( ep, s, t, u ) + 2. * I2D4stu( ep + 1, s, t, u );
      ncomplex sumY = ( 660. * ICache::getI1( ep, Kinem1( msq2 ) ) + 967 * msq2 ) *
                      ( msq2 * msq2 ) * ( msq2 * msq2 ) / 7200.;
      ncomplex sumZ = ( 660. * ICache::getI1( ep, Kinem1( msq1 ) ) + 967 * msq1 ) *
                      ( msq1 * msq1 ) * ( msq1 * msq1 ) / 7200.;
 
      const double ds0tu =
          ( Cay[nss( m, m )] * Cay[nss( n, n )] - Cay[nss( m, n )] * Cay[nss( m, n )] );
      sum1 += sumX * ds0tu;
 
      const double dsvtuY =
          -( Cay[nss( n, n )] - Cay[nss( m, n )] ); /* minus sign of minor v=m */
      sum1 -= sumY * dsvtuY;
 
      const double dsvtuZ =
          +( Cay[nss( m, n )] - Cay[nss( m, m )] ); /* plus sign of minor v=n */
      sum1 -= sumZ * dsvtuZ;
 
      sum1 /= 121 * dstustu;
    }
  }
  else
  {
    assert( ep == 1 );
    const double y11 = Cay[nss( m, m )];
    const double y12 = Cay[nss( m, n )];
    const double y22 = Cay[nss( n, n )];
    sum1 =
        -( y11 * y11 * y11 *
               ( y11 * ( 63 * ( y11 + y12 ) + 7 * y22 ) + 7 * y12 * ( 8 * y12 + 3 * y22 ) +
                 3 * y22 * y22 ) +
           y22 * y22 * y22 *
               ( y22 * ( 63 * ( y22 + y12 ) + 7 * y11 ) + 7 * y12 * ( 8 * y12 + 3 * y11 ) +
                 3 * y11 * y11 ) +
           y12 * y12 * y12 *
               ( 8 * y12 * ( y12 + 3 * y11 + 3 * y22 ) +
                 6 * ( 7 * y11 * y11 + 4 * y11 * y22 + 7 * y22 * y22 ) ) +
           y11 * y22 * y12 * ( 4 * y12 * ( 9 * ( y11 + y22 ) + 4 * y12 ) + 15 * y11 * y22 ) ) /
        2661120;
  }
  pI2D5stu[ep][idx] = sum1;
}
 
/* --------------------------------------------------------
 *   I3D6st triangle in D+12 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I3D6st( int ep, int s, int t ) {
  assert( s != t ); // if (s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I3D6st + ep] ) { I3D6stEval( ep ); }
  int idx = im2( s, t ) - 5;
  return pI3D6st[ep][idx];
}
 
void Minor5::I3D6stEval( int ep ) {
  for ( int s = 1; s <= smax; s++ )
  {
    for ( int t = s + 1; t <= 5; t++ )
    {
      int idx = im2( s, t ) - 5;
 
      const double dstst  = M2( s, t, s, t );
      ncomplex     ivalue = 0;
 
      if ( pmaxS3[idx] <= deps )
      {
        printf( "I3D6 - M2({%d,%d}) <= eps\n", s, t );
        ivalue = std::numeric_limits<double>::quiet_NaN();
      }
      else
      {
        double cf = 1. / 12.;
        for ( int ei = ep; ei <= 1; ei++ )
        {
          ncomplex sum1 = M3( 0, s, t, 0, s, t ) * I3D5st( ei, s, t );
          for ( int u = 1; u <= 5; u++ )
          {
            if ( t == u || s == u ) continue;
            sum1 -= M3( u, s, t, 0, s, t ) * I2D5stu( ei, s, t, u );
          }
          ivalue += cf * sum1;
          cf *= 1. / 6.;
        }
        ivalue /= dstst;
      }
      pI3D6st[ep][idx] = ivalue;
    }
  }
  fEval[E_I3D6st + ep] = true;
}
 
/* --------------------------------------------------------
 *   I2D6stu bubble in D+12 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I2D6stu( int ep, int s, int t, int u ) {
  assert( t != u && u != s && s != t ); // if (t==u || u==s || s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I2D6stu + ep] )
  {
    I2D6stuEval( 0, ep, 1, 2, 3, 4, 5, kinem.p5() );
    I2D6stuEval( 1, ep, 1, 2, 4, 3, 5, kinem.s45() );
    I2D6stuEval( 2, ep, 1, 2, 5, 3, 4, kinem.p4() );
 
    I2D6stuEval( 3, ep, 1, 3, 4, 2, 5, kinem.s12() );
    I2D6stuEval( 4, ep, 1, 3, 5, 2, 4, kinem.s34() );
 
    I2D6stuEval( 5, ep, 1, 4, 5, 2, 3, kinem.p3() );
 
    if ( smax == 5 )
    {
      I2D6stuEval( 6, ep, 2, 3, 4, 1, 5, kinem.p1() );
      I2D6stuEval( 7, ep, 2, 3, 5, 1, 4, kinem.s15() );
 
      I2D6stuEval( 8, ep, 2, 4, 5, 1, 3, kinem.s23() );
 
      I2D6stuEval( 9, ep, 3, 4, 5, 1, 2, kinem.p2() );
    }
 
    fEval[E_I2D6stu + ep] = true;
  }
  int idx = im3( s, t, u ) - 10;
  return pI2D6stu[ep][idx];
}
 
void Minor5::I2D6stuEval( int idx, int ep, int s, int t, int u, int m, int n, double qsq ) {
  ncomplex sum1 = 0;
  if ( ep == 0 )
  {
    const double dstustu  = -2 * qsq; /*M3(s,t,u,s,t,u);*/
    const double msq1     = kinem.mass( m );
    const double msq2     = kinem.mass( n );
    const double s_cutoff = seps1 * pmaxM2[im2( m, n ) - 5];
 
    if ( fabs( dstustu ) <= s_cutoff )
    {
      const double mm12 = msq1 - msq2;
      if ( fabs( mm12 ) < meps )
      {
        sum1 = ( 29 * msq1 + 20. * ICache::getI1( ep, Kinem1( msq1 ) ) ) * ( msq1 * msq1 ) *
               ( msq1 * msq1 ) * msq1 / 14400.;
      }
      else
      {
        const double msq1p6 = ( msq1 * msq1 ) * ( msq1 * msq1 ) * ( msq1 * msq1 );
        const double msq2p6 = ( msq2 * msq2 ) * ( msq2 * msq2 ) * ( msq2 * msq2 );
        sum1                = ( 223 * ( msq1 * msq1p6 - msq2 * msq2p6 ) / 705600 +
                 ( +msq1p6 * ICache::getI1( ep, Kinem1( msq1 ) ) -
                   msq2p6 * ICache::getI1( ep, Kinem1( msq2 ) ) ) /
                     5040. ) /
               mm12;
      }
    }
    else
    {
      ncomplex sumX = 13. * I2D5stu( ep, s, t, u ) + 2. * I2D5stu( ep + 1, s, t, u );
      ncomplex sumY = -( 260. * ICache::getI1( ep, Kinem1( msq2 ) ) + 417 * msq2 ) *
                      ( msq2 * msq2 ) * ( msq2 * msq2 ) * msq2 / 14400.;
      ncomplex sumZ = -( 260. * ICache::getI1( ep, Kinem1( msq1 ) ) + 417 * msq1 ) *
                      ( msq1 * msq1 ) * ( msq1 * msq1 ) * msq1 / 14400.;
 
      const double ds0tu =
          ( Cay[nss( m, m )] * Cay[nss( n, n )] - Cay[nss( m, n )] * Cay[nss( m, n )] );
      sum1 += sumX * ds0tu;
 
      const double dsvtuY =
          -( Cay[nss( n, n )] - Cay[nss( m, n )] ); /* minus sign of minor v=m */
      sum1 -= sumY * dsvtuY;
 
      const double dsvtuZ =
          +( Cay[nss( m, n )] - Cay[nss( m, m )] ); /* plus sign of minor v=n */
      sum1 -= sumZ * dsvtuZ;
 
      sum1 /= 169 * dstustu;
    }
  }
  else
  {
    assert( ep == 1 );
    const double y11 = Cay[nss( m, m )];
    const double y12 = Cay[nss( m, n )];
    const double y22 = Cay[nss( n, n )];
    sum1             = ( y11 * y11 * y11 *
                 ( y11 * ( 21 * y11 * ( 11 * ( y11 + y12 ) + y22 ) + 210 * y12 * y12 +
                           7 * y22 * y22 + 63 * y22 * y12 ) +
                   y12 * y12 * ( 168 * y12 + 112 * y22 ) ) +
             y22 * y22 * y22 *
                 ( y22 * ( 21 * y22 * ( 11 * ( y22 + y12 ) + y11 ) + 210 * y12 * y12 +
                           7 * y11 * y11 + 63 * y11 * y12 ) +
                   y12 * y12 * ( 168 * y12 + 112 * y11 ) ) +
             y12 * y12 *
                 ( y12 * y12 *
                       ( 16 * y12 * y12 + 112 * ( y11 * y11 + y22 * y22 ) +
                         56 * y12 * ( y11 + y22 ) + 120 * y11 * y22 ) +
                   y11 * y22 * ( 90 * y11 * y22 + 140 * y12 * ( y22 + y11 ) ) ) +
             y11 * y11 * y22 * y22 * ( 35 * ( y11 + y22 ) * y12 + 5 * y11 * y22 ) ) /
           138378240;
  }
  pI2D6stu[ep][idx] = sum1;
}
 
/* --------------------------------------------------------
 *   I3D7st triangle in D+14 dim
 * --------------------------------------------------------
 */
ncomplex Minor5::I3D7st( int ep, int s, int t ) {
  assert( s != t ); // if (s==t) return 0;
  if ( ep == 2 ) return 0;
  if ( not fEval[E_I3D7st + ep] ) { I3D7stEval( ep ); }
  int idx = im2( s, t ) - 5;
  return pI3D7st[ep][idx];
}
 
void Minor5::I3D7stEval( int ep ) {
  for ( int s = 1; s <= smax; s++ )
  {
    for ( int t = s + 1; t <= 5; t++ )
    {
      int idx = im2( s, t ) - 5;
 
      const double dstst  = M2( s, t, s, t );
      ncomplex     ivalue = 0;
 
      if ( pmaxS3[idx] <= deps )
      {
        printf( "I3D7 - M2({%d,%d}) <= eps\n", s, t );
        ivalue = std::numeric_limits<double>::quiet_NaN();
      }
      else
      {
        double cf = 1. / 14.;
        for ( int ei = ep; ei <= 1; ei++ )
        {
          ncomplex sum1 = M3( 0, s, t, 0, s, t ) * I3D6st( ei, s, t );
          for ( int u = 1; u <= 5; u++ )
          {
            if ( t == u || s == u ) continue;
            sum1 -= M3( u, s, t, 0, s, t ) * I2D6stu( ei, s, t, u );
          }
          ivalue += cf * sum1;
          cf *= 1. / 7.;
        }
        ivalue /= dstst;
      }
      pI3D7st[ep][idx] = ivalue;
    }
  }
  fEval[E_I3D7st + ep] = true;
}