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-2026 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 The precision of calculation.
27inline constexpr double epsilon = 1e-5;
28//! @brief Alias for the target expression.
29using Expression = std::function<double(const double)>;
30//! @brief Alias for the result of integral.
31using Result = double;
32
33//! @brief Integral methods.
34class Integral
35{
36public:
37 //! @brief Destroy the Integral object.
38 virtual ~Integral() = default;
39
40 //! @brief The operator (()) overloading of Integral class.
41 //! @param lower - lower endpoint
42 //! @param upper - upper endpoint
43 //! @param eps - precision of calculation
44 //! @return result of integral
45 virtual Result operator()(const double lower, const double upper, const double eps) const = 0;
46
47protected:
48 //! @brief Construct a new Integral object.
49 //! @param expr - target expression
50 explicit Integral(Expression expr) : expr{std::move(expr)} {}
51
52 //! @brief Target expression.
53 const Expression expr;
54 //! @brief Calculate the value of the definite integral with the trapezoidal rule.
55 //! @param expr - target expression
56 //! @param left - left endpoint
57 //! @param height - height of trapezoidal
58 //! @param step - number of steps
59 //! @return result of definite integral
60 static double trapezoidalRule(
61 const Expression& expr, const double left, const double height, const std::uint32_t step);
62};
63
64//! @brief The trapezoidal method.
65class Trapezoidal : public Integral
66{
67public:
68 //! @brief Construct a new Trapezoidal object.
69 //! @param expr - target expression
70 explicit Trapezoidal(Expression expr) : Integral(std::move(expr)) {}
71
72 //! @brief The operator (()) overloading of Trapezoidal class.
73 //! @param lower - lower endpoint
74 //! @param upper - upper endpoint
75 //! @param eps - precision of calculation
76 //! @return result of integral
77 Result operator()(const double lower, const double upper, const double eps) const override;
78};
79
80//! @brief The adaptive Simpson's 1/3 method.
81class Simpson : public Integral
82{
83public:
84 //! @brief Construct a new Simpson object.
85 //! @param expr - target expression
86 explicit Simpson(Expression expr) : Integral(std::move(expr)) {}
87
88 //! @brief The operator (()) overloading of Simpson class.
89 //! @param lower - lower endpoint
90 //! @param upper - upper endpoint
91 //! @param eps - precision of calculation
92 //! @return result of integral
93 Result operator()(const double lower, const double upper, const double eps) const override;
94
95private:
96 //! @brief Calculate the value of the definite integral with the Simpson's rule.
97 //! @param left - left endpoint
98 //! @param right - right endpoint
99 //! @param eps - precision of calculation
100 //! @return result of definite integral
101 [[nodiscard]] double simpsonIntegral(const double left, const double right, const double eps) const;
102 //! @brief Composite Simpson's 1/3 formula.
103 //! @param left - left endpoint
104 //! @param right - right endpoint
105 //! @param step - number of steps
106 //! @return result of definite integral
107 [[nodiscard]] double compositeSimpsonOneThird(
108 const double left, const double right, const std::uint32_t step) const;
109 //! @brief Simpson's 1/3 formula.
110 //! @param left - left endpoint
111 //! @param right - right endpoint
112 //! @return result of definite integral
113 [[nodiscard]] double simpsonOneThird(const double left, const double right) const;
114};
115
116//! @brief The Romberg method.
117class Romberg : public Integral
118{
119public:
120 //! @brief Construct a new Romberg object.
121 //! @param expr - target expression
122 explicit Romberg(Expression expr) : Integral(std::move(expr)) {}
123
124 //! @brief The operator (()) overloading of Romberg class.
125 //! @param lower - lower endpoint
126 //! @param upper - upper endpoint
127 //! @param eps - precision of calculation
128 //! @return result of integral
129 Result operator()(const double lower, const double upper, const double eps) const override;
130
131private:
132 //! @brief The Richardson extrapolation.
133 //! @param lowPrec - numerical result obtained with low precision
134 //! @param highPrec - numerical result obtained with high precision
135 //! @param weight - weight factor
136 //! @return extrapolation
137 static double richardsonExtrapolation(const double lowPrec, const double highPrec, const double weight);
138};
139
140//! @brief The Gauss-Legendre's 5-points method.
141class Gauss : public Integral
142{
143public:
144 //! @brief Construct a new Gauss object.
145 //! @param expr - target expression
146 explicit Gauss(Expression expr) : Integral(std::move(expr)) {}
147
148 //! @brief The operator (()) overloading of Gauss class.
149 //! @param lower - lower endpoint
150 //! @param upper - upper endpoint
151 //! @param eps - precision of calculation
152 //! @return result of integral
153 Result operator()(const double lower, const double upper, const double eps) const override;
154
155private:
156 //! @brief Number of Gauss nodes.
157 static constexpr std::uint8_t nodeSize{5};
158 //! @brief Number of Gauss coefficients.
159 static constexpr std::uint8_t coeffSize{2};
160 //! @brief Table of Gauss-Legendre quadrature.
161 static constexpr std::array<std::array<double, coeffSize>, nodeSize> gaussLegendreTbl{
162 ._M_elems: {{-0.9061798459, +0.2369268851},
163 {-0.5384693101, +0.4786286705},
164 {+0.0000000000, +0.5688888889},
165 {+0.5384693101, +0.4786286705},
166 {+0.9061798459, +0.2369268851}}};
167};
168
169//! @brief The Monte-Carlo method.
170class MonteCarlo : public Integral
171{
172public:
173 //! @brief Construct a new MonteCarlo object.
174 //! @param expr - target expression
175 explicit MonteCarlo(Expression expr) : Integral(std::move(expr)) {}
176
177 //! @brief The operator (()) overloading of MonteCarlo class.
178 //! @param lower - lower endpoint
179 //! @param upper - upper endpoint
180 //! @param eps - precision of calculation
181 //! @return result of integral
182 Result operator()(const double lower, const double upper, const double eps) const override;
183
184private:
185 //! @brief Sample from the uniform distribution.
186 //! @param lower - lower endpoint
187 //! @param upper - upper endpoint
188 //! @param eps - precision of calculation
189 //! @return result of definite integral
190 [[nodiscard]] double sampleFromUniformDistribution(const double lower, const double upper, const double eps) const;
191 //! @brief Sample from the normal distribution (Box-Muller transform).
192 //! @param lower - lower endpoint
193 //! @param upper - upper endpoint
194 //! @param eps - precision of calculation
195 //! @return result of definite integral
196 [[deprecated, nodiscard]] double sampleFromNormalDistribution(
197 const double lower, const double upper, const double eps) const;
198};
199} // namespace integral
200} // namespace numeric
201