minoreval.cpp

Go to the documentation of this file.

/*
 * minoreval.cpp - integral coefficient functions
 *
 * this file is part of PJFry library
 * Copyright 2011 Valery Yundin
 */
 
#include "cache.h"
#include "minor.h"
 
/* ========================================================
 * ========================================================
 *
 *                PENTAGON coefficients - 5 point
 *
 * ========================================================
 * ========================================================
 */
 
/* --------------------------------------------------------
    5-point coefficients rank-0
 * --------------------------------------------------------
 */
ncomplex Minor5::evalE( int ep ) {
  ncomplex     ivalue = 0;
  ncomplex     sum1   = 0;
  const double d00    = M1( 0, 0 );
  for ( int s = 1; s <= 5; s++ ) { sum1 += M1( 0, s ) * I4s( ep, s ); }
  ivalue = sum1 / d00;
  return ivalue;
}
 
/* --------------------------------------------------------
    5-point coefficients rank-1
 * --------------------------------------------------------
 */
ncomplex Minor5::evalE( int ep, int i ) {
  ncomplex     ivalue = 0;
  ncomplex     sum1   = 0;
  const double d00    = M1( 0, 0 );
  for ( int s = 1; s <= 5; s++ ) { sum1 += M2( 0, i, 0, s ) * I4s( ep, s ); }
  ivalue = -sum1 / d00;
  return ivalue;
}
 
/* --------------------------------------------------------
    5-point coefficients rank-2
 * --------------------------------------------------------
 */
ncomplex Minor5::evalE( int ep, int i, int j ) {
  ncomplex     ivalue = 0;
  const double d00    = M1( 0, 0 );
  if ( i == 0 && j == 0 )
  {
    ncomplex sum1 = 0;
    for ( int s = 1; s <= 5; s++ ) { sum1 += M1( 0, s ) * I4Ds( ep, s ); }
    ivalue = -sum1 / ( 2 * d00 );
  }
  else
  {
    assert( i != 0 && j != 0 ); // E01, E02, etc do not exist
    ncomplex sum1 = 0;
    for ( int s = 1; s <= 5; s++ )
    {
      sum1 += M2( 0, i, s, j ) * I4Ds( ep, s );
      sum1 += M2( 0, s, 0, j ) * I4Dsi( ep, s, i );
    }
    //     if (i!=j) { // is symmetrization needed?
    //       ncomplex sumX=0;
    //       for (int s=1; s<=5; s++) {
    //         sumX+=M2(0, j, s, i)*I4Ds(ep, s);
    //         sumX+=M2(0, s, 0, i)*I4Dsi(ep, s, j);
    //       }
    //       sum1=0.5*(sum1+sumX);
    //     }
    ivalue = sum1 / d00;
  }
  return ivalue;
}
 
/* --------------------------------------------------------
    5-point coefficients rank-3
 * --------------------------------------------------------
 */
ncomplex Minor5::evalE( int ep, int i, int j, int k ) {
  ncomplex     ivalue = 0;
  const double d00    = M1( 0, 0 );
  if ( i == 0 && j == 0 )
  {
    assert( i == 0 && j == 0 && k != 0 ); // E000 does not exist, E100 is a wrong syntax
 
    /* // Fleischer's formula 6.13
    ncomplex sum1=0;
    for (int s=1; s<=5; s++) {
      sum1+=0.5*M2(0, s, 0, k)*I4Ds(ep, s);
      sum1+=-M1(s, k)*I4D2s(ep, s); // (4-1)/3=1 // NB d-1=3 and d=4
    }
    ivalue=sum1/d00;
    */
    // This variant looks simpler (does not depend on I4D2s)
    ncomplex sum1 = 0;
    for ( int s = 1; s <= 5; s++ )
    {
      sum1 += 0.5 * M2( 0, s, 0, k ) * I4Ds( ep, s );
      sum1 += M1( s, 0 ) * I4D2si( ep, s, k );
    }
    ivalue = sum1 / ( 3 * d00 );
  }
  else
  {
    assert( i != 0 && j != 0 && k != 0 ); // E110, E012, etc do not exist
    ncomplex sum1 = 0;
    ncomplex sumX = 0;
    ncomplex sumY = 0;
    if ( i == j )
    {
      for ( int s = 1; s <= 5; s++ )
      {
        sum1 += 2 * M2( 0, j, s, k ) * I4D2si( ep, s, i );
        sum1 += M2( 0, s, 0, k ) * I4D2sij( ep, s, i, j );
      }
      if ( i != k )
      {
        for ( int s = 1; s <= 5; s++ )
        {
          sumX += M2( 0, j, s, i ) * I4D2si( ep, s, k );
          sumX += M2( 0, k, s, i ) * I4D2si( ep, s, j );
          sumX += M2( 0, s, 0, i ) * I4D2sij( ep, s, k, j );
        }
        sum1 = ( sum1 + 2. * sumX ) / 3.;
      }
    }
    else
    { // i!=j
      for ( int s = 1; s <= 5; s++ )
      {
        sum1 += M2( 0, j, s, k ) * I4D2si( ep, s, i );
        sum1 += M2( 0, i, s, k ) * I4D2si( ep, s, j );
        sum1 += M2( 0, s, 0, k ) * I4D2sij( ep, s, i, j );
      }
      if ( i != k )
      {
        for ( int s = 1; s <= 5; s++ )
        {
          sumX += M2( 0, j, s, i ) * I4D2si( ep, s, k );
          sumX += M2( 0, k, s, i ) * I4D2si( ep, s, j );
          sumX += M2( 0, s, 0, i ) * I4D2sij( ep, s, k, j );
        }
      }
      else { sumX = sum1; }
      if ( j != k )
      {
        for ( int s = 1; s <= 5; s++ )
        {
          sumY += M2( 0, k, s, j ) * I4D2si( ep, s, i );
          sumY += M2( 0, i, s, j ) * I4D2si( ep, s, k );
          sumY += M2( 0, s, 0, j ) * I4D2sij( ep, s, i, k );
        }
      }
      else { sumY = sum1; }
      sum1 = ( sum1 + sumX + sumY ) / 3.;
    }
    ivalue = -sum1 / d00;
  }
  return ivalue;
}
 
/* --------------------------------------------------------
    5-point coefficients rank-4
 * --------------------------------------------------------
 */
ncomplex Minor5::evalE( int ep, int i, int j, int k, int l ) {
  ncomplex     ivalue = 0;
  const double d00    = M1( 0, 0 );
 
  if ( i == 0 && j == 0 )
  {
    if ( k == 0 && l == 0 )
    {
      ncomplex sum1 = 0;
      for ( int s = 1; s <= 5; s++ )
      {
#ifndef USE_GOLEM_MODE
        // Cancel pole and finite part with E00ij - LoopTools-like convention
        sum1 += M1( s, 0 ) * ( I4D2s( ep, s ) + ( ep == 0 ? ( -1. / 6. + 1. / 4. )
                                                          : ( ep == 1 ? -1. / 6. : 0. ) ) );
#else
        // Golem95 convention
        sum1 += M1( s, 0 ) * ( I4D2s( ep, s ) + ( ep == 0 ? -1. / 9. : 0. ) );
#endif
      }
      ivalue = 0.25 * sum1 / d00;
    }
    else
    {
      assert( i == 0 && j == 0 && k != 0 && l != 0 ); // E0001 does not exist, E1200 is a wrong
                                                      // syntax
      ncomplex sum1 = 0;
 
      for ( int s = 1; s <= 5; s++ )
      {
#ifndef USE_GOLEM_MODE
        // Cancel pole and finite part with E0000 - LoopTools-like convention
        sum1 += 0.5 * ( M2( 0, k, s, l ) + M2( 0, l, s, k ) ) *
                ( I4D2s( ep, s ) +
                  ( ep == 0 ? ( -1. / 6. + 1. / 4. ) : ( ep == 1 ? -1. / 6. : 0. ) ) );
#else
        // Golem95 convention
        sum1 += 0.5 * ( M2( 0, k, s, l ) + M2( 0, l, s, k ) ) *
                ( I4D2s( ep, s ) + ( ep == 0 ? -1. / 9. : 0. ) );
#endif
        sum1 += 0.5 * ( M2( 0, s, 0, l ) * I4D2si( ep, s, k ) +
                        M2( 0, s, 0, k ) * I4D2si( ep, s, l ) );
        sum1 += 0.5 * M1( s, 0 ) * ( I4D3sij( ep, s, k, l ) + I4D3sij( ep, s, l, k ) );
      }
      ivalue = -0.25 * sum1 / d00;
    }
  }
  else
  {
    assert( i != 0 && j != 0 && k != 0 && l != 0 ); // E110, E012, etc do not exist
    ncomplex sum1234 = 0;
    for ( int s = 1; s <= 5; s++ )
    {
      ncomplex sum1 = M2( 0, k, s, l ) * I4D3sij( ep, s, i, j ) +
                      M2( 0, i, s, l ) * I4D3sij( ep, s, k, j ) +
                      M2( 0, j, s, l ) * I4D3sij( ep, s, i, k ) +
                      M2( 0, s, 0, l ) * I4D3sijk( ep, s, i, j, k );
      ncomplex sum2 = sum1;
      if ( l != k )
      {
        sum2 = M2( 0, l, s, k ) * I4D3sij( ep, s, i, j ) +
               M2( 0, i, s, k ) * I4D3sij( ep, s, l, j ) +
               M2( 0, j, s, k ) * I4D3sij( ep, s, i, l ) +
               M2( 0, s, 0, k ) * I4D3sijk( ep, s, i, j, l );
      }
      ncomplex sum3 = sum1;
      if ( j == k ) { sum3 = sum2; }
      else if ( l != j )
      {
        sum3 = M2( 0, k, s, j ) * I4D3sij( ep, s, i, l ) +
               M2( 0, i, s, j ) * I4D3sij( ep, s, k, l ) +
               M2( 0, l, s, j ) * I4D3sij( ep, s, i, k ) +
               M2( 0, s, 0, j ) * I4D3sijk( ep, s, i, l, k );
      }
      ncomplex sum4 = sum1;
      if ( i == j ) { sum4 = sum3; }
      else if ( l != i )
      {
        sum4 = M2( 0, k, s, i ) * I4D3sij( ep, s, l, j ) +
               M2( 0, l, s, i ) * I4D3sij( ep, s, k, j ) +
               M2( 0, j, s, i ) * I4D3sij( ep, s, l, k ) +
               M2( 0, s, 0, i ) * I4D3sijk( ep, s, l, j, k );
      }
      sum1234 += sum1 + sum2 + sum3 + sum4;
    }
    ivalue = sum1234 / ( 4 * d00 );
  }
  return ivalue;
}
 
/* --------------------------------------------------------
    5-point coefficients rank-5
 * --------------------------------------------------------
 */
ncomplex Minor5::evalE( int ep, int i, int j, int k, int l, int m ) {
  ncomplex     ivalue = 0;
  const double d00    = M1( 0, 0 );
 
  if ( i == 0 && j == 0 )
  {
    if ( k == 0 && l == 0 )
    {
      assert( m != 0 ); // E00000 does not exist
      ncomplex sum1 = 0;
      for ( int s = 1; s <= 5; s++ )
      {
        sum1 += +M2( 0, s, 0, m ) * I4D2s( ep, s ) +
                M1( s, 0 ) * ( 4. * I4D3si( ep, s, m ) -
                               2. * ( ep < 2 ? I4D3si( ep + 1, s, m ) : 0. ) );
      }
      ivalue = -sum1 / ( 20. * d00 );
    }
    else
    {
      assert( i == 0 && j == 0 && k != 0 && l != 0 && m != 0 ); // E00012 does not exist,
                                                                // E00100 is a wrong syntax
      ncomplex sum1 = 0;
      ncomplex sumX = 0;
      ncomplex sumY = 0;
      if ( k == l )
      {
        for ( int s = 1; s <= 5; s++ )
        {
          sum1 += +2 * M2( 0, k, s, m ) * I4D3si( ep, s, l ) +
                  M2( 0, s, 0, m ) * I4D3sij( ep, s, k, l );
        }
        if ( ep == 0 ) sum1 += 1. / 24. * ( M1( k, m ) - M2( 0, k, l, m ) );
        if ( k != m )
        {
          for ( int s = 1; s <= 5; s++ )
          {
            sumX += +M2( 0, m, s, k ) * I4D3si( ep, s, l ) +
                    M2( 0, l, s, k ) * I4D3si( ep, s, m ) +
                    M2( 0, s, 0, k ) * I4D3sij( ep, s, m, l );
          }
          if ( ep == 0 ) sumX += 1. / 48. * ( M1( m, k ) + M1( l, k ) - M2( 0, l, m, k ) );
          sum1 = ( sum1 + 2. * sumX ) / 3.;
        }
      }
      else
      { // k!=l
        for ( int s = 1; s <= 5; s++ )
        {
          sum1 += +M2( 0, k, s, m ) * I4D3si( ep, s, l ) +
                  M2( 0, l, s, m ) * I4D3si( ep, s, k ) +
                  M2( 0, s, 0, m ) * I4D3sij( ep, s, k, l );
        }
        if ( ep == 0 )
          sum1 += 1. / 48. * ( M1( k, m ) + M1( l, m ) - M2( 0, k, l, m ) - M2( 0, l, k, m ) );
        if ( k != m )
        {
          for ( int s = 1; s <= 5; s++ )
          {
            sumX += +M2( 0, m, s, k ) * I4D3si( ep, s, l ) +
                    M2( 0, l, s, k ) * I4D3si( ep, s, m ) +
                    M2( 0, s, 0, k ) * I4D3sij( ep, s, m, l );
          }
          if ( ep == 0 )
            sumX +=
                1. / 48. * ( M1( m, k ) + M1( l, k ) - M2( 0, m, l, k ) - M2( 0, l, m, k ) );
        }
        else { sumX = sum1; }
        if ( l != m )
        {
          for ( int s = 1; s <= 5; s++ )
          {
            sumY += +M2( 0, k, s, l ) * I4D3si( ep, s, m ) +
                    M2( 0, m, s, l ) * I4D3si( ep, s, k ) +
                    M2( 0, s, 0, l ) * I4D3sij( ep, s, k, m );
          }
          if ( ep == 0 )
            sumY +=
                1. / 48. * ( M1( k, l ) + M1( m, l ) - M2( 0, k, m, l ) - M2( 0, m, k, l ) );
        }
        else { sumY = sum1; }
        sum1 = ( sum1 + sumX + sumY ) / 3.;
      }
      sumX = 0;
      for ( int s = 1; s <= 5; s++ ) { sumX += M1( s, 0 ) * I4D4sijk( ep, s, k, l, m ); }
      sum1   = 3. * sum1 + 2. * sumX;
      ivalue = sum1 / ( 10. * d00 );
    }
  }
  else
  {
    assert( i != 0 && j != 0 && k != 0 && l != 0 && m != 0 );
    ncomplex sum12345 = 0;
    for ( int s = 1; s <= 5; s++ )
    {
      ncomplex sum1 = +M2( 0, l, s, m ) * I4D4sijk( ep, s, i, j, k ) +
                      M2( 0, i, s, m ) * I4D4sijk( ep, s, l, j, k ) +
                      M2( 0, j, s, m ) * I4D4sijk( ep, s, i, l, k ) +
                      M2( 0, k, s, m ) * I4D4sijk( ep, s, i, j, l ) +
                      M2( 0, s, 0, m ) * I4D4sijkl( ep, s, i, j, k, l );
      ncomplex sum2 = sum1;
      if ( m != l )
      {
        sum2 = +M2( 0, m, s, l ) * I4D4sijk( ep, s, i, j, k ) +
               M2( 0, i, s, l ) * I4D4sijk( ep, s, m, j, k ) +
               M2( 0, j, s, l ) * I4D4sijk( ep, s, i, m, k ) +
               M2( 0, k, s, l ) * I4D4sijk( ep, s, i, j, m ) +
               M2( 0, s, 0, l ) * I4D4sijkl( ep, s, i, j, k, m );
      }
      ncomplex sum3 = sum1;
      if ( k == l ) { sum3 = sum2; }
      else if ( m != k )
      {
        sum3 = +M2( 0, l, s, k ) * I4D4sijk( ep, s, i, j, m ) +
               M2( 0, i, s, k ) * I4D4sijk( ep, s, l, j, m ) +
               M2( 0, j, s, k ) * I4D4sijk( ep, s, i, l, m ) +
               M2( 0, m, s, k ) * I4D4sijk( ep, s, i, j, l ) +
               M2( 0, s, 0, k ) * I4D4sijkl( ep, s, i, j, m, l );
      }
      ncomplex sum4 = sum1;
      if ( j == k ) { sum4 = sum3; }
      else if ( m != j )
      {
        sum4 = +M2( 0, l, s, j ) * I4D4sijk( ep, s, i, m, k ) +
               M2( 0, i, s, j ) * I4D4sijk( ep, s, l, m, k ) +
               M2( 0, m, s, j ) * I4D4sijk( ep, s, i, l, k ) +
               M2( 0, k, s, j ) * I4D4sijk( ep, s, i, m, l ) +
               M2( 0, s, 0, j ) * I4D4sijkl( ep, s, i, m, k, l );
      }
      ncomplex sum5 = sum1;
      if ( i == j ) { sum5 = sum4; }
      else if ( m != i )
      {
        sum5 = +M2( 0, l, s, i ) * I4D4sijk( ep, s, m, j, k ) +
               M2( 0, m, s, i ) * I4D4sijk( ep, s, l, j, k ) +
               M2( 0, j, s, i ) * I4D4sijk( ep, s, m, l, k ) +
               M2( 0, k, s, i ) * I4D4sijk( ep, s, m, j, l ) +
               M2( 0, s, 0, i ) * I4D4sijkl( ep, s, m, j, k, l );
      }
      sum12345 += sum1 + sum2 + sum3 + sum4 + sum5;
    }
    ivalue = -sum12345 / ( 5 * d00 );
  }
  return ivalue;
}
 
/* ========================================================
 * ========================================================
 *
 *                BOX coefficients - 4 point
 *
 * ========================================================
 * ========================================================
 */
#ifdef USE_GOLEM_MODE_6
int psix;
#  define ix( i ) i < pm5->psix ? i : i - 1
#else
#  define ix( i ) i
#endif
/* --------------------------------------------------------
    4-point scalar for GolemMode
 * --------------------------------------------------------
 */
#ifdef USE_GOLEM_MODE
ncomplex Minor4::A( int ep ) {
  ncomplex ivalue = pm5->I4s( ep, ps );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
/* --------------------------------------------------------
    4-point coefficients rank-1
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor4::A( int ep, int i ) {
  ncomplex ivalue = -pm5->I4Dsi( ep, ps, ix( i + offs ) );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor4::evalD( int ep, int i ) {
  ncomplex ivalue = 0;
  if ( i >= ps || ps == 5 ) i = i + 1;
  ivalue = -pm5->I4Dsi( ep, ps, i );
  return ivalue;
}
 
/* --------------------------------------------------------
    4-point coefficients rank-2
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor4::A( int ep, int i, int j ) {
  ncomplex ivalue = pm5->I4D2sij( ep, ps, ix( i + offs ), ix( j + offs ) );
  return ivalue;
}
 
ncomplex Minor4::B( int ep ) {
  ncomplex ivalue = -0.5 * pm5->I4Ds( ep, ps );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor4::evalD( int ep, int i, int j ) {
  ncomplex ivalue = 0;
 
  if ( i == 0 && j == 0 ) { ivalue = -0.5 * pm5->I4Ds( ep, ps ); }
  else
  {
    assert( i != 0 && j != 0 ); // D01, D02, etc do not exist
    if ( i >= ps || ps == 5 ) i = i + 1;
    if ( j >= ps || ps == 5 ) j = j + 1;
 
    ivalue = pm5->I4D2sij( ep, ps, i, j );
  }
  return ivalue;
}
 
/* --------------------------------------------------------
    4-point coefficients rank-3
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor4::A( int ep, int i, int j, int k ) {
  ncomplex ivalue = -pm5->I4D3sijk( ep, ps, ix( i + offs ), ix( j + offs ), ix( k + offs ) );
  return ivalue;
}
 
ncomplex Minor4::B( int ep, int k ) {
  ncomplex ivalue = 0.5 * pm5->I4D2si( ep, ps, ix( k + offs ) );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor4::evalD( int ep, int i, int j, int k ) {
  ncomplex ivalue = 0;
 
  if ( i == 0 && j == 0 )
  {
    assert( k != 0 ); // D000 does not exist, D100 is a wrong syntax
    if ( k >= ps || ps == 5 ) k = k + 1;
 
    ivalue = 0.5 * pm5->I4D2si( ep, ps, k );
  }
  else
  {
    assert( i != 0 && j != 0 && k != 0 ); // D110, D012, etc do not exist
    if ( i >= ps || ps == 5 ) i = i + 1;
    if ( j >= ps || ps == 5 ) j = j + 1;
    if ( k >= ps || ps == 5 ) k = k + 1;
    ivalue = -pm5->I4D3sijk( ep, ps, i, j, k );
  }
  return ivalue;
}
 
/* --------------------------------------------------------
    4-point coefficients rank-4
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor4::A( int ep, int i, int j, int k, int l ) {
  ncomplex ivalue =
      pm5->I4D4sijkl( ep, ps, ix( i + offs ), ix( j + offs ), ix( k + offs ), ix( l + offs ) );
  return ivalue;
}
 
ncomplex Minor4::B( int ep, int k, int l ) {
  ncomplex ivalue = -0.5 * pm5->I4D3sij( ep, ps, ix( k + offs ), ix( l + offs ) );
  return ivalue;
}
 
ncomplex Minor4::C( int ep ) {
  ncomplex ivalue = 0.25 * pm5->I4D2s( ep, ps );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor4::evalD( int ep, int i, int j, int k, int l ) {
  ncomplex ivalue = 0;
  if ( i == 0 && j == 0 )
  {
    if ( k == 0 && l == 0 ) { ivalue = 0.25 * pm5->I4D2s( ep, ps ); }
    else
    {
      assert( i == 0 && j == 0 && k != 0 && l != 0 ); // D0001 does not exist, D1200 is a wrong
                                                      // syntax
      if ( k >= ps || ps == 5 ) k = k + 1;
      if ( l >= ps || ps == 5 ) l = l + 1;
 
      ivalue = -0.5 * pm5->I4D3sij( ep, ps, k, l );
    }
  }
  else
  {
    assert( i != 0 && j != 0 && k != 0 && l != 0 ); // D110, D012, etc do not exist
    if ( i >= ps || ps == 5 ) i = i + 1;
    if ( j >= ps || ps == 5 ) j = j + 1;
    if ( k >= ps || ps == 5 ) k = k + 1;
    if ( l >= ps || ps == 5 ) l = l + 1;
    ivalue = pm5->I4D4sijkl( ep, ps, i, j, k, l );
  }
  return ivalue;
}
 
/* ========================================================
 * ========================================================
 *
 *                TRIANGLE coefficients - 3 point
 *
 * ========================================================
 * ========================================================
 */
/* --------------------------------------------------------
    3-point scalar for GolemMode
 * --------------------------------------------------------
 */
#ifdef USE_GOLEM_MODE
ncomplex Minor3::A( int ep ) {
  ncomplex ivalue = pm5->I3st( ep, ps, pt );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
/* --------------------------------------------------------
    3-point coefficients rank-1
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor3::A( int ep, int i ) {
  ncomplex ivalue = -pm5->I3Dsti( ep, ps, pt, ix( i + offs ) );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor3::evalC( int ep, int i ) {
  ncomplex ivalue = 0;
  if ( i >= ps ) i = i + 1;
  if ( i >= pt || pt == 5 )
  {
    i = i + 1;
    if ( i == ps ) i = i + 1;
  }
  ivalue = -pm5->I3Dsti( ep, ps, pt, i );
  return ivalue;
}
 
/* --------------------------------------------------------
    3-point coefficients rank-2
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor3::A( int ep, int i, int j ) {
  ncomplex ivalue = pm5->I3D2stij( ep, ps, pt, ix( i + offs ), ix( j + offs ) );
  return ivalue;
}
 
ncomplex Minor3::B( int ep ) {
  ncomplex ivalue = -0.5 * pm5->I3Dst( ep, ps, pt );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor3::evalC( int ep, int i, int j ) {
  ncomplex ivalue = 0;
 
  if ( i == 0 && j == 0 ) { ivalue = -0.5 * pm5->I3Dst( ep, ps, pt ); }
  else
  {
    assert( i != 0 && j != 0 ); // C01, C02, etc do not exist
    int tmp;
    tswap( i, j, tmp );
 
    if ( i >= ps ) i = i + 1;
    if ( i >= pt || pt == 5 )
    {
      i = i + 1;
      if ( i == ps ) i = i + 1;
    }
    if ( j >= ps ) j = j + 1;
    if ( j >= pt || pt == 5 )
    {
      j = j + 1;
      if ( j == ps ) j = j + 1;
    }
 
    ivalue = pm5->I3D2stij( ep, ps, pt, i, j );
  }
  return ivalue;
}
 
/* --------------------------------------------------------
    3-point coefficients rank-3
 * --------------------------------------------------------
 */
// Global chord indexing, golem-like
#ifdef USE_GOLEM_MODE
ncomplex Minor3::A( int ep, int i, int j, int k ) {
  ncomplex ivalue =
      -pm5->I3D3stijk( ep, ps, pt, ix( i + offs ), ix( j + offs ), ix( k + offs ) );
  return ivalue;
}
 
ncomplex Minor3::B( int ep, int k ) {
  ncomplex ivalue = 0.5 * pm5->I3D2sti( ep, ps, pt, ix( k + offs ) );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor3::evalC( int ep, int i, int j, int k ) {
  ncomplex ivalue = 0;
 
  if ( i == 0 && j == 0 )
  {
    assert( i == 0 && j == 0 && k != 0 ); // C000 does not exist, C100 is a wrong syntax
    if ( k >= ps ) k = k + 1;
    if ( k >= pt || pt == 5 )
    {
      k = k + 1;
      if ( k == ps ) k = k + 1;
    }
 
    ivalue = 0.5 * pm5->I3D2sti( ep, ps, pt, k );
  }
  else
  {
    assert( i != 0 && j != 0 && k != 0 ); // D110, D012, etc do not exist
    if ( i >= ps ) i = i + 1;
    if ( i >= pt || pt == 5 )
    {
      i = i + 1;
      if ( i == ps ) i = i + 1;
    }
    if ( j >= ps ) j = j + 1;
    if ( j >= pt || pt == 5 )
    {
      j = j + 1;
      if ( j == ps ) j = j + 1;
    }
    if ( k >= ps ) k = k + 1;
    if ( k >= pt || pt == 5 )
    {
      k = k + 1;
      if ( k == ps ) k = k + 1;
    }
    ivalue = -pm5->I3D3stijk( ep, ps, pt, i, j, k );
  }
  return ivalue;
}
 
/* ========================================================
 * ========================================================
 *
 *             BUBBLE coefficients - 2 point
 *
 * ========================================================
 * ========================================================
 */
/* --------------------------------------------------------
    2-point scalar for GolemMode
 * --------------------------------------------------------
 */
#ifdef USE_GOLEM_MODE
ncomplex Minor2::A( int ep ) {
  ncomplex ivalue = pm5->I2stu( ep, ps, pt, pu );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
/* --------------------------------------------------------
    2-point coefficients rank-1
 * --------------------------------------------------------
 */
#ifdef USE_GOLEM_MODE
ncomplex Minor2::A( int ep, int i ) {
  ncomplex ivalue = -pm5->I2Dstui( ep, ps, pt, pu, ix( i + offs ) );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor2::evalB( int ep, int i ) {
  ncomplex ivalue = 0;
 
  if ( i >= ps ) i = i + 1;
  if ( i >= pt )
  {
    i = i + 1;
    if ( i == ps ) i = i + 1;
  }
  if ( i >= pu || pu == 5 )
  {
    i = i + 1;
    if ( i == ps ) i = i + 1;
    if ( i == pt ) i = i + 1;
  }
  ivalue = -pm5->I2Dstui( ep, ps, pt, pu, i );
  return ivalue;
}
 
/* --------------------------------------------------------
    2-point coefficients rank-2
 * --------------------------------------------------------
 */
#ifdef USE_GOLEM_MODE
ncomplex Minor2::A( int ep, int i, int j ) {
  ncomplex ivalue = pm5->I2D2stuij( ep, ps, pt, pu, ix( i + offs ), ix( j + offs ) );
  return ivalue;
}
 
ncomplex Minor2::B( int ep ) {
  ncomplex ivalue = -0.5 * pm5->I2Dstu( ep, ps, pt, pu );
  return ivalue;
}
#endif /* USE_GOLEM_MODE */
 
// Local chord indexing
ncomplex Minor2::evalB( int ep, int i, int j ) {
  ncomplex ivalue = 0;
 
  if ( i == 0 && j == 0 ) { ivalue = -0.5 * pm5->I2Dstu( ep, ps, pt, pu ); }
  else
  {
    assert( i != 0 && j != 0 ); // B01, B02, etc do not exist
    if ( i >= ps ) i = i + 1;
    if ( i >= pt )
    {
      i = i + 1;
      if ( i == ps ) i = i + 1;
    }
    if ( i >= pu || pu == 5 )
    {
      i = i + 1;
      if ( i == ps ) i = i + 1;
      if ( i == pt ) i = i + 1;
    }
    ivalue = pm5->I2D2stuij( ep, ps, pt, pu, i, i );
  }
  return ivalue;
}