61 numIterations_(numIterations),
62 weight_(numIterations + 1),
63 root_(numIterations + 1)
65 calculateWeightAndRoot();
74 const std::vector<double>&
getWeight() const noexcept
85 const std::vector<double>&
getRoot() const noexcept
98 double derivative{ 0.0 };
107 Result(
const double val,
const double deriv) noexcept :
117 void calculateWeightAndRoot() noexcept
119 const double numIterationsDouble =
static_cast<double>(numIterations_);
120 for (
uint32 step = 0; step <= numIterations_; ++step)
122 double root =
std::cos(
constants::pi * (
static_cast<double>(step) - 0.25) / (numIterationsDouble + 0.5));
123 Result result = calculatePolynomialValueAndDerivative(root);
125 double newtonRaphsonRatio;
128 newtonRaphsonRatio = result.value / result.derivative;
129 root -= newtonRaphsonRatio;
130 result = calculatePolynomialValueAndDerivative(root);
131 }
while (std::fabs(newtonRaphsonRatio) > EPSILON);
134 weight_[step] = 2.0 / ((1.0 -
utils::sqr(root)) * result.derivative * result.derivative);
145 Result calculatePolynomialValueAndDerivative(
const double x) noexcept
147 Result result(x, 0.0);
149 double value_minus_1 = 1.0;
151 for (
uint32 step = 2; step <= numIterations_; ++step)
153 const double stepDouble =
static_cast<double>(step);
154 const double value = ((2.0 * stepDouble - 1.0) * x * result.value - (stepDouble - 1.0) * value_minus_1) / stepDouble;
155 result.derivative = stepDouble *
f * (x * value - result.value);
157 value_minus_1 = result.value;
158 result.value = value;
165 const double EPSILON{ 1
e-15 };
167 const uint32 numIterations_;
168 std::vector<double> weight_;
169 std::vector<double> root_;
184 const std::function<
double(
double)>&
f)
187 const std::vector<double>& weight = legendrePolynomial.
getWeight();
188 const std::vector<double>& root = legendrePolynomial.
getRoot();
190 const double width = 0.5 * (high - low);
191 const double mean = 0.5 * (low + high);
193 double gaussLegendre = 0.0;
194 for (
uint32 step = 1; step <= n; ++step)
196 gaussLegendre += weight[step] *
f(width * root[step] +
mean);
199 return gaussLegendre * width;