#ifndef EIGEN_MONEGATO_H
#define EIGEN_MONEGATO_H
namespace Eigen
{
/**
* \ingroup NumericalIntegration_Module
*
* \class Monegato
*
* \brief A class for computing the Gauss/Kronrod weights.
*
* \tparam Scalar floating point type
*
* This class is based on "Some remarks on the construction of
* extended Gaussian quadrature rules", Giovanni Monegato,
* Math. Comp., Vol. 32 (1978) pp. 247-252. http://www.jstor.org/stable/2006272 .
*
* The code is based on quadpackcpp (http://quadpackpp.sourceforge.net/)
* library written by Jerry Gagelman <jerry@os-scientific.org> which is distributed under
* GNU General Public License
*/
template <typename Scalar>
class Monegato
{
public:
/**
* \brief Compute the absolute value of a scalar
* \param[in] x inpute scalar value
* \return absolute value of scalar input
*/
static Scalar abs(Scalar x)
{
return (x < Scalar(0)) ? -(x) : x;
}
/**
* \brief Compute the Legendre polynomials and their error
* \param[in] n degree
* \param[in] x value of the variable
* \param[out] err error to be returned
* \return value of the Legendre polynomial
*
* This function is based on the routine gsl_sf_legendre_Pl_e
* distributed with GSL
*/
static Scalar legendre_err(const int n, const Scalar x, Scalar& err)
{
if (n == 0)
{
err = Scalar(0);
return Scalar(1);
}
else if (n == 1)
{
err = Scalar(0);
return x;
}
// below Sree modified this to avoid -Wmaybe-uninitialized
Scalar P0 = Scalar(1), P1 = x, P2=x;
Scalar E0 = NumTraits<Scalar>::epsilon();
Scalar E1 = abs(x) * NumTraits<Scalar>::epsilon();
for (int k = 1; k < n; ++k)
{
P2 = ((2*k + 1) * x * P1 - k * P0) / (k + 1);
err = ((2*k + 1) * abs(x) * E1 + k * E0) / (2*(k + 1));
P0 = P1; P1 = P2;
E0 = E1; E1 = err;
}
return P2;
}
/**
* \brief Three-term recursion identity for the Legendre derivatives
* \param[in] n degree
* \param[in] x value of the variable
* \return value of the Legendre derivative
*/
static Scalar legendre_deriv(int const n, Scalar const x)
{
if (n == 0)
{
return Scalar(0);
}
else if (n == 1)
{
return Scalar(1);
}
Scalar P0 = Scalar(1);
Scalar P1 = x;
Scalar P2 = P1;
Scalar dP0 = Scalar(0);
Scalar dP1 = Scalar(1);
Scalar dP2 = P1;
for (int k = 1; k < n; ++k)
{
P2 = ((2*k + 1) * x * P1 - k * P0) / (k + Scalar(1));
dP2 = (2*k + 1) * P1 + dP0;
P0 = P1;
P1 = P2;
dP0 = dP1;
dP1 = dP2;
}
return dP2;
}
/**
* \brief Compute derivatives of Chebyshev polynomials
* \param[in] x variable
* \param[in] n_ degree
* \param[in] coefs Chebyshev coefficients
* \return Derivative of Chebyshev polynomials
*/
static Scalar chebyshev_series_deriv(const Scalar x,
const int n_,
const Eigen::Array<Scalar,
Eigen::Dynamic, 1>& coefs)
{
Scalar d1(0), d2(0);
Scalar y2 = 2 * x; // linear term for Clenshaw recursion
for (int k = n_; k >= 2; --k)
{
Scalar temp = d1;
d1 = y2 * d1 - d2 + k * coefs(k);
d2 = temp;
}
return y2 * d1 - d2 + coefs(1);
}
/**
* \breif Evaluation of the Chebyshev polynomial using Clenshaw recursion
* \param[in] x variable
* \param[in] n_ degree
* \param[in] coefs Chebyshev coefficients
* \param[out] err error to be returned
* \return value of Chebyshev polynomial
*/
static Scalar chebyshev_series(const Scalar x,
const int n_,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& coefs,
Scalar& err)
{
Scalar d1(0), d2(0);
Scalar absc = abs(coefs(0)); // final term for truncation error
Scalar y2 = 2 * x; // linear term for Clenshaw recursion
for (int k = n_; k >= 1; --k)
{
Scalar temp = d1;
d1 = y2 * d1 - d2 + coefs(k);
d2 = temp;
absc += abs(coefs(k));
}
err = absc * NumTraits<Scalar>::epsilon();
return x * d1 - d2 + coefs(0)/2.;
}
/**
* \brief Computes the zeros of the Legendre polynomial
* \param[in] m_ degree
* \param[out] zeros the zeros of the Legendre polynomial
*/
static void legendre_zeros(const int m_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& zeros)
{
Eigen::Array<Scalar, Eigen::Dynamic, 1> temp = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(m_+1);
zeros(0) = Scalar(-1);
zeros(1) = Scalar(1);
Scalar delta, epsilon;
for (int k = 1; k <= m_; ++k)
{
// Loop to locate zeros of P_k interlacing z_0,...,z_k
for (int j = 0; j < k; ++j)
{
// Newton's method for P_k :
// initialize solver at midpoint of (z_j, z_{j+1})
delta = 1;
Scalar x_j = (zeros(j) + zeros(j+1)) / 2.;
Scalar P_k = legendre_err(k, x_j, epsilon);
while (abs(P_k) > epsilon &&
abs(delta) > NumTraits<Scalar>::epsilon())
{
delta = P_k / legendre_deriv(k, x_j);
x_j -= delta;
P_k = legendre_err(k, x_j, epsilon);
}
temp(j) = x_j;
}
// Copy roots tmp_0,...,tmp_{k-1} to z_1,...z_k:
zeros(k+1) = zeros(k);
for (int j = 0; j < k; ++j)
{
zeros(j+1) = temp(j);
}
}
}
/**
* \brief Computes coefficients of the Chebyshev polynomial.
* \param[in] m_ degree of the Chebyshev polynomial.
* \param[out] coefs coefficients of the Chebyshev polynomial.
*/
static void chebyshev_coefs(const int m_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& coefs)
{
int ell = (m_ + 1)/2;
Eigen::Array<Scalar, Eigen::Dynamic, 1> alpha = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(ell+1);
Eigen::Array<Scalar, Eigen::Dynamic, 1> f = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(ell+1);
// Care must be exercised in initalizing the constants in the definitions.
// Compilers interpret expressions like "(2*k + 1.0)/(k + 2.0)" as floating
// point precision, before casting to Real.
f(1) = Scalar(m_+1) / Scalar(2*m_ + 3);
alpha(0) = Scalar(1); // coefficient of T_{m+1}
alpha(1) = -f(1);
for (int k = 2; k <= ell; ++k)
{
f(k) = f(k-1) * (2*k - 1) * (m_ + k) / (k * (2*m_ + 2*k + 1));
alpha(k) = -f(k);
for (int i = 1; i < k; ++i)
{
alpha(k) -= f(i) * alpha(k-i);
}
}
for (int k = 0; k <= ell; ++k)
{
coefs(m_ + 1 - 2*k) = alpha(k);
if (m_ >= 2*k)
{
coefs(m_ - 2*k) = Scalar(0);
}
}
}
/**
* \brief Compute the Gauss/Kronrod abscissae
* \param[in] n_ size of Gauss-Kronrod arrays
* \param[in] Gauss-Legendre degree
* \param[in] coefs Chebyshev coefficients
* \param[out] xgk_ Gauss/Kronrod abscissae
*/
static void gauss_kronrod_abscissae(const int n_,
const int m_,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& zeros,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& coefs,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& xgk_)
{
// Now from the function gauss_kronrod_abscissae
Scalar delta, epsilon;
for (int k = 0; k < n_ / 2; ++k)
{
delta = 1;
// Newton's method for E_{n+1} :
Scalar x_k = (zeros(m_-k) + zeros(m_+1-k))/Scalar(2);
Scalar E = chebyshev_series(x_k,n_,coefs, epsilon);
while (abs(E) > epsilon &&
abs(delta) > NumTraits<Scalar>::epsilon() )
{
delta = E / chebyshev_series_deriv(x_k,n_,coefs);
x_k -= delta;
E = chebyshev_series(x_k,n_,coefs, epsilon);
}
xgk_(2*k) = x_k;
// Copy adjacent Legendre-zero into the array:
if (2*k+1 < n_)
{
xgk_(2*k+1) = zeros(m_-k);
}
}
}
/**
* \brief Compute Gauss/Kronrod weights
* \param[in] n_ size of Gauss-Kronrod arrays
* \param[in] m_ Gauss-Legendre degree
* \param[in] xgk_ Gauss/Kronrod abscissae
* \param[out] wg_ Gauss weights
* \param[out] wgk_ Kronrod weights
*/
static void gauss_kronrod_weights(const int& n_,
const int m_,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& coefs,
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& xgk_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& wg_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& wgk_)
{
Scalar err;
// Gauss-Legendre weights:
for(int k = 0; k < n_ / 2; ++k)
{
Scalar x = xgk_(2*k + 1);
wg_(k) = (Scalar(-2) / ((m_ + 1) * legendre_deriv(m_, x) * legendre_err(m_+1, x, err)));
}
// The ratio of leading coefficients of P_n and T_{n+1} is computed
// from the recursive formulae for the respective polynomials.
Scalar F_m = Scalar(2) / Scalar(2*m_ + 1);
for (int k = 1; k <= m_; ++k)
{
F_m *= (Scalar(2*k) / Scalar(2*k - 1));
}
// Gauss-Kronrod weights:
for (int k = 0; k < n_; ++k)
{
Scalar x = xgk_(k);
if (k % 2 == 0)
{
wgk_(k) = F_m / (legendre_err(m_, x, err) * chebyshev_series_deriv(x,n_,coefs));
}
else
{
wgk_(k) = (wg_(k/2) + F_m / (legendre_deriv(m_, x) * chebyshev_series(x,n_,coefs, err)));
}
}
}
/**
* \brief A method ofr computing Gauss/Kronrod abscissae and weights
* \param[in] m_ Gauss-Legendre degree
* \param[out] xgk_ Gauss/Kronrod abscissae
* \param[out] wgk_ Gauss/Kronrod weights
* \param[out] xk_ Gauss abscissae
* \param[out] wg_ Gauss weights
*
* Note that Gauss abscissae is not calculated here.
*/
static void computeAbscissaeAndWeights(unsigned int m_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& xgk_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& wgk_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& xk_,
Eigen::Array<Scalar, Eigen::Dynamic, 1>& wg_)
{
const unsigned int n_ = m_ + 1;
xgk_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_); //2*nNodes+1
wgk_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_); //2*nNodes+1
xk_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_/2); //2*nNodes
wg_ = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_/2); //2*nNodes
// initialise the coefficients to zero
Eigen::Array<Scalar, Eigen::Dynamic, 1> coefs = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n_+1);
Eigen::Array<Scalar, Eigen::Dynamic, 1> zeros = Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(m_+2);
legendre_zeros(m_, zeros);
chebyshev_coefs(m_, coefs);
gauss_kronrod_abscissae(n_, m_, zeros,coefs, xgk_);
gauss_kronrod_weights(n_, m_,coefs, xgk_, wg_, wgk_ );
}
};
} // namespace Eigen
#endif // EIGEN_MONEGATO_H