Page MenuHomeHEPForge
Files F9501254

GaussLegendreQuadratureQ.icc
No OneTemporary
Actions

GaussLegendreQuadratureQ.icc
View Options

#include <cmath>
#include <stdexcept>
#include <algorithm>
#include "npstat/nm/PreciseType.hh"
namespace npstat {
template <typename FcnResult, typename FcnArg>
inline lapack_double GaussLegendreQuadratureQ::integrate(
const Functor1<FcnResult,FcnArg>& function,
const lapack_double left, const lapack_double right) const
{
typedef PreciseType<double>::type precise_type;
if (buf_.size() != npoints_)
buf_.resize(npoints_);
std::pair<lapack_double, lapack_double>* results = &buf_[0];
const unsigned halfpoints = npoints_/2;
const lapack_double midpoint = (left + right)/2.0;
const lapack_double unit = (right - left)/2.0;
FcnArg a;
lapack_double v;
for (unsigned i=0; i<halfpoints; ++i)
{
a = static_cast<FcnArg>(midpoint + unit*a_[i]);
v = w_[i]*static_cast<lapack_double>(function(a));
results[2*i].first = std::abs(v);
results[2*i].second = v;
a = static_cast<FcnArg>(midpoint - unit*a_[i]);
v = w_[i]*static_cast<lapack_double>(function(a));
results[2*i+1].first = std::abs(v);
results[2*i+1].second = v;
}
std::sort(buf_.begin(), buf_.end());
precise_type sum = 0.0;
for (unsigned i=0; i<npoints_; ++i)
sum += results[i].second;
return sum*unit;
}
template <typename FcnResult, typename FcnArg>
inline lapack_double GaussLegendreQuadratureQ::integrate(
const Functor1<FcnResult,FcnArg>& function,
const lapack_double left, const lapack_double right,
const unsigned nsplit) const
{
typedef PreciseType<double>::type precise_type;
if (!nsplit) throw std::invalid_argument(
"In npstat::GaussLegendreQuadratureQ::integrate: "
"number of subintervals must be positive");
const lapack_double step = (right - left)/nsplit;
precise_type acc = 0.0;
lapack_double b = left;
for (unsigned i=0; i<nsplit; ++i)
{
const lapack_double a = b;
b = (i == nsplit - 1U ? right : a + step);
acc += integrate(function, a, b);
}
return acc;
}
template <class Functor>
inline std::vector<std::pair<lapack_double,lapack_double> >
GaussLegendreQuadratureQ::weightedIntegrationPoints(
const Functor& weight,
const lapack_double left, const lapack_double right) const
{
typedef std::pair<lapack_double,lapack_double> DDPair;
std::vector<DDPair> result;
result.reserve(npoints_);
const unsigned halfpoints = npoints_/2;
const lapack_double midpoint = (left + right)/2.0;
const lapack_double unit = (right - left)/2.0;
for (unsigned i=0; i<halfpoints; ++i)
{
lapack_double x = midpoint + unit*a_[i];
lapack_double w = unit*w_[i]*weight(x);
result.push_back(DDPair(x, w));
x = midpoint - unit*a_[i];
w = unit*w_[i]*weight(x);
result.push_back(DDPair(x, w));
}
return result;
}
}

File Metadata

Mime Type
text/x-c
Expires
Sun, Feb 23, 2:09 PM (2 h, 18 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486475
Default Alt Text
GaussLegendreQuadratureQ.icc (3 KB)

Event Timeline