1//! @file integral.hpp
2//! @author ryftchen
3//! @brief The declarations (integral) in the numeric module.
4//! @version 0.1.0
5//! @copyright Copyright (c) 2022-2025 ryftchen. All rights reserved.
6
7#pragma once
8
9#include <cstdint>
10#include <functional>
11
12//! @brief The numeric module.
13namespace numeric // NOLINT(modernize-concat-nested-namespaces)
14{
15//! @brief Integral-related functions in the numeric module.
16namespace integral
17{
18//! @brief Brief function description.
19//! @return function description (module_function)
20inline static const char* description() noexcept
21{
22 return "NUM_INTEGRAL";
23}
24extern const char* version() noexcept;
25
26//! @brief Target expressions.
27using Expression = std::function<double(const double)>;
28//! @brief The precision of calculation.
29inline constexpr double epsilon = 1e-5;
30
31//! @brief Integral methods.
32class Integral
33{
34public:
35 //! @brief Construct a new Integral object.
36 Integral() = default;
37 //! @brief Destroy the Integral object.
38 virtual ~Integral() = default;
39 //! @brief Construct a new Integral object.
40 Integral(const Integral&) = default;
41 //! @brief Construct a new Integral object.
42 Integral(Integral&&) noexcept = default;
43 //! @brief The operator (=) overloading of Integral class.
44 //! @return reference of the Integral object
45 Integral& operator=(const Integral&) = default;
46 //! @brief The operator (=) overloading of Integral class.
47 //! @return reference of the Integral object
48 Integral& operator=(Integral&&) noexcept = default;
49
50 //! @brief The operator (()) overloading of Integral class.
51 //! @param lower - lower endpoint
52 //! @param upper - upper endpoint
53 //! @param eps - precision of calculation
54 //! @return result of integral
55 virtual double operator()(const double lower, const double upper, const double eps) const = 0;
56
57protected:
58 //! @brief Calculate the value of the definite integral with the trapezoidal rule.
59 //! @param expr - target expression
60 //! @param left - left endpoint
61 //! @param height - height of trapezoidal
62 //! @param step - number of steps
63 //! @return result of definite integral
64 static double trapezoidalRule(
65 const Expression& expr, const double left, const double height, const std::uint32_t step);
66};
67
68//! @brief The trapezoidal method.
69class Trapezoidal : public Integral
70{
71public:
72 //! @brief Construct a new Trapezoidal object.
73 //! @param expr - target expression
74 explicit Trapezoidal(Expression expr) : expr{std::move(expr)} {}
75
76 //! @brief The operator (()) overloading of Trapezoidal class.
77 //! @param lower - lower endpoint
78 //! @param upper - upper endpoint
79 //! @param eps - precision of calculation
80 //! @return result of integral
81 double operator()(const double lower, const double upper, const double eps) const override;
82
83private:
84 //! @brief Target expression.
85 const Expression expr;
86};
87
88//! @brief The adaptive Simpson's 1/3 method.
89class Simpson : public Integral
90{
91public:
92 //! @brief Construct a new Simpson object.
93 //! @param expr - target expression
94 explicit Simpson(Expression expr) : expr{std::move(expr)} {}
95
96 //! @brief The operator (()) overloading of Simpson class.
97 //! @param lower - lower endpoint
98 //! @param upper - upper endpoint
99 //! @param eps - precision of calculation
100 //! @return result of integral
101 double operator()(const double lower, const double upper, const double eps) const override;
102
103private:
104 //! @brief Target expression.
105 const Expression expr;
106 //! @brief Calculate the value of the definite integral with the Simpson's rule.
107 //! @param left - left endpoint
108 //! @param right - right endpoint
109 //! @param eps - precision of calculation
110 //! @return result of definite integral
111 [[nodiscard]] double simpsonIntegral(const double left, const double right, const double eps) const;
112 //! @brief Composite Simpson's 1/3 formula.
113 //! @param left - left endpoint
114 //! @param right - right endpoint
115 //! @param step - number of steps
116 //! @return result of definite integral
117 [[nodiscard]] double compositeSimpsonOneThird(
118 const double left, const double right, const std::uint32_t step) const;
119 //! @brief Simpson's 1/3 formula.
120 //! @param left - left endpoint
121 //! @param right - right endpoint
122 //! @return result of definite integral
123 [[nodiscard]] double simpsonOneThird(const double left, const double right) const;
124};
125
126//! @brief The Romberg method.
127class Romberg : public Integral
128{
129public:
130 //! @brief Construct a new Romberg object.
131 //! @param expr - target expression
132 explicit Romberg(Expression expr) : expr{std::move(expr)} {}
133
134 //! @brief The operator (()) overloading of Romberg class.
135 //! @param lower - lower endpoint
136 //! @param upper - upper endpoint
137 //! @param eps - precision of calculation
138 //! @return result of integral
139 double operator()(const double lower, const double upper, const double eps) const override;
140
141private:
142 //! @brief Target expression.
143 const Expression expr;
144 //! @brief The Richardson extrapolation.
145 //! @param lowPrec - numerical result obtained with low precision
146 //! @param highPrec - numerical result obtained with high precision
147 //! @param weight - weight factor
148 //! @return extrapolation
149 static double richardsonExtrapolation(const double lowPrec, const double highPrec, const double weight);
150};
151
152//! @brief The Gauss-Legendre's 5-points method.
153class Gauss : public Integral
154{
155public:
156 //! @brief Construct a new Gauss object.
157 //! @param expr - target expression
158 explicit Gauss(Expression expr) : expr{std::move(expr)} {}
159
160 //! @brief The operator (()) overloading of Gauss class.
161 //! @param lower - lower endpoint
162 //! @param upper - upper endpoint
163 //! @param eps - precision of calculation
164 //! @return result of integral
165 double operator()(const double lower, const double upper, const double eps) const override;
166
167private:
168 //! @brief Target expression.
169 const Expression expr;
170 //! @brief Number of Gauss nodes.
171 static constexpr std::uint8_t nodeSize{5};
172 //! @brief Number of Gauss coefficients.
173 static constexpr std::uint8_t coeffSize{2};
174 //! @brief Table of Gauss-Legendre quadrature.
175 static constexpr std::array<std::array<double, coeffSize>, nodeSize> gaussLegendreTbl{
176 ._M_elems: {{-0.9061798459, +0.2369268851},
177 {-0.5384693101, +0.4786286705},
178 {+0.0000000000, +0.5688888889},
179 {+0.5384693101, +0.4786286705},
180 {+0.9061798459, +0.2369268851}}};
181};
182
183//! @brief The Monte-Carlo method.
184class MonteCarlo : public Integral
185{
186public:
187 //! @brief Construct a new MonteCarlo object.
188 //! @param expr - target expression
189 explicit MonteCarlo(Expression expr) : expr{std::move(expr)} {}
190
191 //! @brief The operator (()) overloading of MonteCarlo class.
192 //! @param lower - lower endpoint
193 //! @param upper - upper endpoint
194 //! @param eps - precision of calculation
195 //! @return result of integral
196 double operator()(const double lower, const double upper, const double eps) const override;
197
198private:
199 //! @brief Target expression.
200 const Expression expr;
201 //! @brief Sample from the uniform distribution.
202 //! @param lower - lower endpoint
203 //! @param upper - upper endpoint
204 //! @param eps - precision of calculation
205 //! @return result of definite integral
206 [[nodiscard]] double sampleFromUniformDistribution(const double lower, const double upper, const double eps) const;
207 //! @brief Sample from the normal distribution (Box-Muller transform).
208 //! @param lower - lower endpoint
209 //! @param upper - upper endpoint
210 //! @param eps - precision of calculation
211 //! @return result of definite integral
212 [[deprecated, nodiscard]] double sampleFromNormalDistribution(
213 const double lower, const double upper, const double eps) const;
214};
215} // namespace integral
216} // namespace numeric
217