1//! @file integral.cpp
2//! @author ryftchen
3//! @brief The definitions (integral) in the numeric module.
4//! @version 0.1.0
5//! @copyright Copyright (c) 2022-2025 ryftchen. All rights reserved.
6
7#include "integral.hpp"
8
9#include <random>
10
11namespace numeric::integral
12{
13//! @brief Function version number.
14//! @return version number (major.minor.patch)
15const char* version() noexcept
16{
17 static const char* const ver = "0.1.0";
18 return ver;
19}
20
21double Integral::trapezoidalRule(
22 const Expression& expr, const double left, const double height, const std::uint32_t step)
23{
24 double sum = 0.0;
25 double x = left;
26 const double delta = height / step;
27 for (std::uint32_t i = 0; i < step; ++i)
28 {
29 const double area = (expr(x) + expr(x + delta)) * delta / 2.0;
30 sum += area;
31 x += delta;
32 }
33 return sum;
34}
35
36double Trapezoidal::operator()(const double lower, const double upper, const double eps) const
37{
38 if (!expr)
39 {
40 return 0.0;
41 }
42
43 const double height = upper - lower;
44 double s1 = 0.0;
45 double s2 = 0.0;
46 std::uint32_t n = 1;
47 do
48 {
49 const double sum = trapezoidalRule(expr, left: lower, height, step: n);
50 s1 = s2;
51 s2 = sum;
52 n *= 2;
53 }
54 while (std::fabs(x: s1 - s2) > eps);
55 return s2;
56}
57
58double Simpson::operator()(const double lower, const double upper, const double eps) const
59{
60 return expr ? simpsonIntegral(left: lower, right: upper, eps) : 0.0;
61}
62
63double Simpson::simpsonIntegral(const double left, const double right, const double eps) const
64{
65 const double mid = (left + right) / 2.0;
66 const double sum = simpsonOneThird(left, right);
67 if (std::fabs(x: sum - (compositeSimpsonOneThird(left, right: mid, step: 2) + compositeSimpsonOneThird(left: mid, right, step: 2))) > eps)
68 {
69 return simpsonIntegral(left, right: mid, eps) + simpsonIntegral(left: mid, right, eps);
70 }
71 return sum;
72}
73
74double Simpson::compositeSimpsonOneThird(const double left, const double right, const std::uint32_t step) const
75{
76 const double stepLen = (right - left) / step;
77 double sum = 0.0;
78 for (std::uint32_t i = 0; i < step; ++i)
79 {
80 sum += simpsonOneThird(left: left + (i * stepLen), right: left + ((i + 1) * stepLen));
81 }
82 return sum;
83}
84
85double Simpson::simpsonOneThird(const double left, const double right) const
86{
87 return (expr(left) + 4.0 * expr((left + right) / 2.0) + expr(right)) / 6.0 * (right - left);
88}
89
90double Romberg::operator()(const double lower, const double upper, const double eps) const
91{
92 if (!expr)
93 {
94 return 0.0;
95 }
96
97 const double height = upper - lower;
98 const auto trapezoid = std::bind(f&: trapezoidalRule, args: std::ref(t: expr), args: lower, args: height, args: std::placeholders::_1);
99 std::uint32_t k = 0;
100
101 std::vector<std::vector<double>> rombergTbl{};
102 rombergTbl.emplace_back(args: std::vector<double>{trapezoid(1)});
103 do
104 {
105 ++k;
106 rombergTbl.emplace_back();
107 rombergTbl[k].reserve(n: k + 1);
108 rombergTbl[k].emplace_back(args: trapezoid(std::pow(x: 2, y: k)));
109
110 for (std::uint32_t i = 1; i <= k; ++i)
111 {
112 rombergTbl[k].emplace_back(
113 args: richardsonExtrapolation(lowPrec: rombergTbl[k - 1][i - 1], highPrec: rombergTbl[k][i - 1], weight: std::pow(x: 4, y: i)));
114 }
115 }
116 while (std::fabs(x: rombergTbl[k][k] - rombergTbl[k - 1][k - 1]) > eps);
117 return rombergTbl[k][k];
118}
119
120double Romberg::richardsonExtrapolation(const double lowPrec, const double highPrec, const double weight)
121{
122 return (weight * highPrec - lowPrec) / (weight - 1.0);
123}
124
125double Gauss::operator()(const double lower, const double upper, const double eps) const
126{
127 if (!expr)
128 {
129 return 0.0;
130 }
131
132 double s1 = 0.0;
133 double s2 = 0.0;
134 std::uint32_t n = 1;
135 do
136 {
137 double sum = 0.0;
138 const double stepLen = (upper - lower) / n;
139 for (std::uint32_t i = 0; i < n; ++i)
140 {
141 const double left = lower + (i * stepLen);
142 const double right = left + stepLen;
143 for (const auto& coeff : gaussLegendreTbl)
144 {
145 const double x = ((right - left) * coeff[0] + (left + right)) / 2.0;
146 const double polynomial = expr(x) * coeff[1] * (right - left) / 2.0;
147 sum += polynomial;
148 }
149 }
150 s1 = s2;
151 s2 = sum;
152 n *= 2;
153 }
154 while (std::fabs(x: s1 - s2) > eps);
155 return s2;
156}
157
158double MonteCarlo::operator()(const double lower, const double upper, const double eps) const
159{
160 return expr ? sampleFromUniformDistribution(lower, upper, eps) : 0.0;
161}
162
163double MonteCarlo::sampleFromUniformDistribution(const double lower, const double upper, const double eps) const
164{
165 const std::uint32_t n = 1.0 / eps;
166 std::mt19937_64 engine(std::random_device{}());
167 std::uniform_real_distribution<double> dist(lower, upper);
168 double sum = 0.0;
169 for (std::uint32_t i = 0; i < n; ++i)
170 {
171 sum += expr(dist(engine));
172 }
173 sum *= (upper - lower) / n;
174 return sum;
175}
176
177double MonteCarlo::sampleFromNormalDistribution(const double lower, const double upper, const double eps) const
178{
179 const std::uint32_t n = 1.0 / eps;
180 const double mu = (lower + upper) / 2.0;
181 const double sigma = (upper - lower) / 6.0;
182 std::mt19937_64 engine(std::random_device{}());
183 std::uniform_real_distribution<double> dist(0.0, 1.0);
184 double sum = 0.0;
185 double x = 0.0;
186 for (std::uint32_t i = 0; i < n; ++i)
187 {
188 do
189 {
190 const double u1 = dist(engine);
191 const double u2 = dist(engine);
192 const double mag = sigma * std::sqrt(x: -2.0 * std::log(x: u1));
193 x = mag * std::sin(x: 2.0 * std::numbers::pi * u2) + mu;
194 }
195 while ((x < lower) || (x > upper));
196 const double probDens = (1.0 / std::sqrt(x: 2.0 * std::numbers::pi * sigma * sigma))
197 * std::pow(x: std::numbers::e, y: (-(x - mu) * (x - mu)) / (2.0 * sigma * sigma));
198 sum += expr(x) / probDens;
199 }
200 sum /= n;
201 return sum;
202}
203} // namespace numeric::integral
204