#ifndef EIGEN_PIESSENS_H
#define EIGEN_PIESSENS_H
namespace Eigen
{
/**
* \ingroup NumericalIntegration_Module
*
* \class Piessens
*
* \brief This class computes Kronrod abscissae & weights for arbitrary precision
*
* \tparam Scalar floating point type
*
* This class is based on the work by R. Piessens, et.al,published in the
* journal Mathematics of Computation, Volume 28, Number 125, January, 1974.
*/
template <typename Scalar>
class Piessens
{
public:
/**
* \brief kronrod adds n+1 points to an n-point Gaussian rule.
*
* This function is a C++ implementation of original work by R. Piessens, et.al,
* published in the journal Mathematics of Computation, Volume 28, Number 125,
* January, 1974. Where possible strucutre and variable naming convention
* has been aligned to favor work of the QUADPACK Gaus-Kronrod integration
* routines created by individuals of the same group.
*
* This function calculates the abscissas and weights of the (2n+1)-point
* Gauss-Kronrod quadrature formula which is obtained from the n-point
* Gauss quadrature formula with the optimal addition of (n+1)-points.
*
* The optimally added points are the Kronrod abscissae. The
* abscissas and weights for both the Gauss and Gauss Kronrod rules
* are calculated for integration over the interval (-1, +1).
*
* Because the quadrature formula is symmetric with respect to the origin,
* only the positive abscissas are calculated. Weights corresponding to the
* symetric abscissae are equal. Weights of weightGauss are calculated as well.
*
* Work by Dr. John Burkhardt made note that the code published in Mathematics of
* Computation omitted the definition of the second Chebyshev coefficient (chebCoeff2),
* and Dr. Burkhardt's contributions are reflected here with permission.
*
* The arrays abcsGaussKronrod, weightGaussKronrod and weightGauss contain the
* positive abscissae in decreasing order, and the weights of each abscissa in
* the Gauss-Kronrod and Gauss rules, respectively.
*
* Ported to C++/Eigen and templated for multiprecision by Mark Sauder,
* Sreekumar Thaithara Balan, Matt Beall, and R. Jeff Jenkins - September 2014.
*
* Input Parameters:
* \param[in] n, the order of the Gauss rule.
*
* Return Parameters:
* \param[in,out] abscGaussKronrod[n+1] The Gauss-Kronrod abscissae.
* \param[in,out] weightGaussKronrod[n+1] The weights for the Gauss-Kronrod rule.
* \param[in,out] weightGauss[n+1] The weights for the Gauss rule.
*/
static void kronrod(unsigned int nNodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& abscGaussKronrod,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGaussKronrod,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& weightGauss)
{
unsigned int arraySize = nNodes + 1;
abscGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);
weightGaussKronrod = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize);
weightGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(arraySize / 2);
Scalar aN(0.0);
Scalar d(2.0);
for (size_t i = 0; i < nNodes; ++i)
{
aN += Scalar(1.0);
d *= aN / (aN + Scalar(0.5));
}
unsigned int m = (nNodes + 1) / 2;
bool even = (nNodes == 2 * m);
// aK is an index variable to account for calculating only the positive abscissae
Scalar aK = aN;
// Calculation of the Chebyshev coefficients of the orthogonal polynomial.
Eigen::Array<Scalar, Eigen::Dynamic, 1> tau(m);
tau(0) = (aN + Scalar(2.0)) / (Scalar(2) * aN + Scalar(3.0));
Eigen::Array<Scalar, Eigen::Dynamic, 1> betaCoeffs(m + 1);
betaCoeffs(m - 1) = tau(0) - Scalar(1.0);
for (size_t k = 1; k < m; ++k)
{
// This step accounts for both positive and negative abscissae
aK += Scalar(2.0);
tau(k) = ((aK - Scalar(1.0)) * aK - aN * (aN + Scalar(1.0) )) * (aK + Scalar(2.0) ) * tau(k-1) /
(aK * ((aK + Scalar(3.0) ) * (aK + Scalar(2.0) ) - aN * (aN + Scalar(1.0) )));
betaCoeffs(m-k-1) = tau(k);
for (size_t i = 1; i <= k; ++i)
{
betaCoeffs(m - k - 1) = betaCoeffs(m - k + i - 1) * tau(i - 1) + betaCoeffs(m - k - 1);
}
}
betaCoeffs(m) = Scalar(1.);
// Calculation of approximate values for the abscissae as inital values
// for the Newton-Raphson iterative solution. These values are derived
// from Pythagorean identities to the original code to more closely follow
// the mathematics of the 1974 CoM paper.
// @TODO The usage of constant Pi with fixed precision needs to be changed to the following for multiprecision
//RealScalar pi = NumTraits<RealScalar>::Pi();
Scalar s1 = sin((M_PI / Scalar(2) ) / (Scalar(2.) * aN + Scalar(1.0) ));
Scalar c1 = cos((M_PI / Scalar(2) ) / (Scalar(2.) * aN + Scalar(1.0) ));
Scalar s2 = sin((M_PI) / (Scalar(2.) * aN + Scalar(1.0) ));
Scalar c2 = cos((M_PI) / (Scalar(2.) * aN + Scalar(1.0) ));
// Coefficient for Gauss and Kronrod abscissae and weights
Scalar chebCoeff1 = Scalar(1.0) - Scalar(1.0) / (Scalar(8.0) * aN * aN) + Scalar(1.0) / (Scalar(8.0) * aN * aN * aN);
Scalar chebCoeff2 = Scalar(2.0) / (Scalar(2. * nNodes + 1));
for (size_t i = 1; i <= nNodes; ++i)
{
chebCoeff2 = Scalar(4.0) * chebCoeff2 * i / (nNodes + i);
}
Scalar abscK = chebCoeff1 * c1;
Scalar temp(0.);
// Calculation of the K-th (Kronrod) abscissa and the corresponding weight.
for (size_t k = 0; k < nNodes; ++k)
{
abscWeightKronrod(nNodes, m, even, chebCoeff2, betaCoeffs, abscK, weightGaussKronrod(k));
abscGaussKronrod(k) = abscK;
++k;
temp = c1;
c1 = temp * c2 - s1 * s2;
s1 = temp * s2 + s1 * c2;
abscK = chebCoeff1 * c1;
// Calculation of the k+1 (Gauss) abscissa and the corresponding weights.
abscWeightGauss(nNodes, m, even, chebCoeff2, betaCoeffs, abscK, weightGaussKronrod(k),
weightGauss(k/2));
abscGaussKronrod(k) = abscK;
temp = c1;
c1 = temp * c2 - s1 * s2;
s1 = temp * s2 + s1 * c2;
abscK = chebCoeff1 * c1;
}
// Add a Kronrod abscissa at the origin if n is even.
if (even)
{
abscWeightKronrod(nNodes, m, even, chebCoeff2, betaCoeffs, abscK, weightGaussKronrod(nNodes));
}
// Set the abscissa value at the origin to zero and exit the function.
abscGaussKronrod(nNodes) = Scalar(0.0);
return;
}
/**
* \brief abscWeightKronrod calculates a Kronrod abscissa and weight.
*
* Input Parameters:
* \param[in] betaCoeffs[m+1] The Chebyshev coefficients.
* \param[in] chebCoeff A value needed to compute weights.
* \param[in] even A boolean variable that is TRUE if n is even.
* \param[in] n The order of the Gauss rule.
* \param[in] m The value of ( n + 1 ) / 2.
*
* Input/output:
* \param[in,out] abscGaussKronrod An estimate for the abscissa on input and the computed abscissa on output.
* \param[in,out] weightGaussKronrod The Gauss-Kronrod weight.
*/
static void abscWeightKronrod(unsigned int nNodes,
unsigned int m,
bool even,
Scalar chebCoeff,
Eigen::Array<Scalar,
Eigen::Dynamic, 1> betaCoeffs,
Scalar& abscGaussKronrod,
Scalar& weightGaussKronrod)
{
Scalar ai;
Scalar b0(0.);
Scalar b1(0.);
Scalar b2(0.);
Scalar d0(0.);
Scalar d1(0.);
Scalar d2(0.);
Scalar delta(1.);
Scalar dif(0.);
Scalar f(0.);
Scalar fd(0.);
Scalar yy(0.);
int i = 0;
size_t iter = 0;
size_t iterationLimit = 50;
// Iterative process for the computation of a Kronrod abscissa.
while (abs(delta) > machineEpsilon())
{
++iter;
b1 = Scalar(0.0);
b2 = betaCoeffs(m);
yy = Scalar(4.0) * abscGaussKronrod * abscGaussKronrod - Scalar(2.0);
d1 = Scalar(0.0);
if (even)
{
ai = Scalar(m + m + 1);
d2 = ai * betaCoeffs(m);
dif = Scalar(2.0);
}
else
{
ai = Scalar(m + 1);
d2 = Scalar(0.0);
dif = Scalar(1.0);
}
for (size_t k = 0; k < m; ++k)
{
ai -= dif;
i = m - (int)k - 1;
b0 = b1;
b1 = b2;
d0 = d1;
d1 = d2;
b2 = yy * b1 - b0 + betaCoeffs(i);
if (!even)
{
i += 1;
}
d2 = yy * d1 - d0 + ai * betaCoeffs(i);
}
if (even)
{
f = abscGaussKronrod * (b2 - b1);
fd = d2 + d1;
}
else
{
f = Scalar(0.5) * (b2 - b0);
fd = Scalar(4.0) * abscGaussKronrod * d2;
}
// Newton correction.
delta = f / fd;
abscGaussKronrod -= delta;
if (abscGaussKronrod == Scalar(0.))
{
break;
}
// Identify non-convergence of the iterative solver after 50 iterations
if (iter > iterationLimit)
{
std::cout << "@abscWeightKronrod Newton-Raphson iterative abscissae solver failed."<<std::endl;
return;
}
}
// Computation of the weight.
d0 = Scalar(1.);
d1 = abscGaussKronrod;
ai = Scalar(0.);
for (size_t k = 0; k < nNodes - 1; ++k)
{
ai = ai + Scalar(1.);
d2 = ((ai + ai + Scalar(1.)) * (abscGaussKronrod * d1) - (ai * d0)) / (ai + Scalar(1.));
d0 = d1;
d1 = d2;
}
weightGaussKronrod = chebCoeff / (fd * d2);
return;
}
/**
* \brief abscWeightGauss calculates a Gaussian abscissa and two weights.
*
* Input Parameters
* \param[in] betaCoeffs[m+1] The Chebyshev coefficients.
* \param[in] chebCoeff A value needed to compute weights.
* \param[in] even A boolean variable that is TRUE if n is even.
* \param[in] nNodes The order of the Gauss rule.
* \param[in] m The value of ( nNodes + 1 ) / 2.
*
* Input/Output
* \param[in,out] abscGaussKronrod An estimate for the abscissa on input and the computed abscissa on output.
* \param[in,out] weightGaussKronrod The Gauss-Kronrod weight.
* \param[in,out] weightGauss The Gauss weight.
*/
static void abscWeightGauss(unsigned int nNodes,
unsigned int m,
bool even,
const Scalar& chebCoeff,
Eigen::Array<Scalar, Eigen::Dynamic, 1> betaCoeffs,
Scalar& abscGaussKronrod,
Scalar& weightGaussKronrod,
Scalar& weightGauss)
{
Scalar ai(0.);
Scalar delta(1.);
Scalar p0(0.);
Scalar p1(0.);
Scalar p2(0.);
Scalar pd0(0.);
Scalar pd1(0.);
Scalar pd2(0.);
Scalar yy(0.);
size_t iter = 0;
size_t iterationLimit = 50;
// Iterative process for the computation of a Gaussian abscissa.
while (abs(delta) > machineEpsilon())
{
++iter;
p0 = Scalar(1.);
p1 = abscGaussKronrod;
pd0 = Scalar(0.);
pd1 = Scalar(1.);
// If nNodes <= 1, initialize p2 and pd2 to avoid problems calculating delta.
if (nNodes <= 1)
{
if (machineEpsilon() < abs(abscGaussKronrod))
{
p2 = (Scalar(3.0) * (abscGaussKronrod) * (abscGaussKronrod) - Scalar(1.0)) / Scalar(2.0);
pd2 = Scalar(3.0) * (abscGaussKronrod);
}
else
{
p2 = Scalar(3.0) * (abscGaussKronrod);
pd2 = Scalar(3.0);
}
}
ai = Scalar(0.0);
for (size_t k = 0; k < nNodes - 1; ++k)
{
ai = ai + Scalar(1.0);
p2 = ((ai + ai + Scalar(1.0) ) * abscGaussKronrod * p1 - ai * p0) / (ai + Scalar(1.0) );
pd2 = ((ai + ai + Scalar(1.0) ) * (p1 + abscGaussKronrod * pd1) - ai * pd0) / (ai + Scalar(1.0) );
p0 = p1;
p1 = p2;
pd0 = pd1;
pd1 = pd2;
}
// Newton correction.
delta = p2 / pd2;
if (abscGaussKronrod == Scalar(0.0))
{
abscGaussKronrod -= delta;
break;
}
else
{
abscGaussKronrod -= delta;
}
// Identify non-convergence of the iterative solver after iteration limit.
if (iter > iterationLimit)
{
std::cout << "@abscWeightGauss Newton-Raphson iterative abscissae solver failed."<<std::endl;;
return;
}
}
// Computation of the Gauss weight.
Scalar aN = nNodes;
weightGauss = Scalar(2.0) / (aN * pd2 * p0);
// Initialize
p1 = Scalar(0.0);
p2 = betaCoeffs(m);
yy = Scalar(4.0) * (abscGaussKronrod) * (abscGaussKronrod) - Scalar(2.0);
for (size_t k = 1; k <= m; ++k)
{
p0 = p1;
p1 = p2;
p2 = yy * p1 - p0 + betaCoeffs(m - k);
}
if (even)
{
weightGaussKronrod = weightGauss + chebCoeff / (abscGaussKronrod * pd2 * (p2 - p1));
}
else
{
weightGaussKronrod = weightGauss + (Scalar(2.0) * chebCoeff) / (pd2 * (p2 - p0));
}
return;
}
/**
* \brief machineEpsilon returns the machine precision roundoff error.
* \returns Returns the machine precision round-off error for double precision
*/
static Scalar machineEpsilon()
{
// ISO standard machine precision values
// half precision: 2^-10 (9.76563e-04)
// single precision: 2^-23 (1.1920928955078125e-07)
// double precision: 2^-52 (2.22044604925031308e-16)
// extended precision: 2^-63 (1.08420217248550443e-19)
// quad precision: 2^-112 (1.9259299e-34)
return NumTraits<Scalar>::epsilon();
}
/**
* \brief A method ofr computing Gauss/Kronrod abscissae and weights
* \param nNodes Gauss-Legendre degree
* \param abscGaussKronrod Gauss/Kronrod abscissae
* \param weightGaussKronrod Gauss/Kronrod weights
* \param abscGauss Gauss abscissae
* \param weightGauss Gauss weights
*/
static void computeAbscissaeAndWeights(unsigned int nNodes,
Eigen::Array<Scalar, Eigen::Dynamic, 1> &abscGaussKronrod,
Eigen::Array<Scalar, Eigen::Dynamic, 1> &weightGaussKronrod,
Eigen::Array<Scalar, Eigen::Dynamic, 1> abscGauss,
Eigen::Array<Scalar, Eigen::Dynamic, 1> &weightGauss)
{
Piessens::kronrod(nNodes,abscGaussKronrod,weightGaussKronrod,weightGauss);
abscGauss = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(abscGauss.rows());
}
};
}//namespace Eigen
#endif //EIGEN_PIESSENS_H