adaptive_simpson.hpp¶
按误差估计自适应细分区间，数值计算平滑函数的定积分；达到深度上限或浮点分辨率时停止，误差估计不保证实际误差上界。
Complexity: Time: \(O(2^{d})\) function evaluations in the worst case, where \(d\) is lim. Space: \(O(d)\) recursion stack.
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(2^d) function evaluations in the worst case, where d is lim.
/// Space: O(d) recursion stack.
#include <cassert>
#include <cmath>
#include <functional>
#include <numeric>
#include <utility>
namespace noya {
/// @brief Numerically integrate a smooth real function on [l, r] with
/// adaptive Simpson subdivision. eps targets the estimated absolute error;
/// recursion also stops at lim or floating-point resolution, without a certified bound.
template <class Function, class Real = long double>
Real adaptive_simpson(Function fun, Real l, Real r, Real eps = Real(1e-10),
int lim = 24) {
assert(eps > Real{});
assert(lim >= 0);
if (l == r) {
return Real{};
}
bool rev = r < l;
if (rev) std::swap(l, r);
// Keep the usual arithmetic when finite; split before adding on overflow.
auto midpoint = [](Real a, Real b) {
Real s = a + b;
return std::isfinite(s) ? s / Real(2) : a / Real(2) + b / Real(2);
};
auto sim = [](Real a, Real b, Real fa, Real fm, Real fb) {
Real wid = b - a;
Real val = wid * (fa + Real(4) * fm + fb) / Real(6);
if (std::isfinite(val)) return val;
Real avg = std::lerp(std::midpoint(fa, fb), fm, Real(2) / Real(3));
return std::isfinite(wid) ? wid * avg
: ((b / Real(2) - a / Real(2)) * avg) * Real(2);
};
Real mdd = midpoint(l, r);
Real fl = std::invoke(fun, l);
Real fm = std::invoke(fun, mdd);
Real fr = std::invoke(fun, r);
Real sum = sim(l, r, fl, fm, fr);
auto sol = [&](auto &dfs, Real a, Real b, Real mid, Real fa, Real fc,
Real fb, Real est, Real tol, int dep) -> Real {
if (mid == a || mid == b) return (b - a) * std::midpoint(fa, fb);
Real lm = midpoint(a, mid);
Real rm = midpoint(mid, b);
Real flm = std::invoke(fun, lm);
Real frm = std::invoke(fun, rm);
Real sl = sim(a, mid, fa, flm, fc);
Real sr = sim(mid, b, fc, frm, fb);
Real dif = sl + sr - est;
if (dep == 0 || std::abs(dif) <= Real(15) * tol) {
return sl + sr + dif / Real(15);
}
return dfs(dfs, a, mid, lm, fa, flm, fc, sl, tol / Real(2), dep - 1) +
dfs(dfs, mid, b, rm, fc, frm, fb, sr, tol / Real(2), dep - 1);
};
Real res = sol(sol, l, r, mdd, fl, fm, fr, sum, eps, lim);
return rev ? -res : res;
}
} // namespace noya
// NOYA_ADAPTIVE_SIMPSON_HPP
#ifndef NOYA_ADAPTIVE_SIMPSON_HPP
#define NOYA_ADAPTIVE_SIMPSON_HPP 1
/// @complexity Time: O(2^d) function evaluations in the worst case, where d is lim.
/// Space: O(d) recursion stack.
#include <cassert>
#include <cmath>
#include <functional>
#include <numeric>
#include <utility>
namespace noya {
/// @brief Numerically integrate a smooth real function on [l, r] with
/// adaptive Simpson subdivision. eps targets the estimated absolute error;
/// recursion also stops at lim or floating-point resolution, without a certified bound.
template <class Function, class Real = long double>
Real adaptive_simpson(Function fun, Real l, Real r, Real eps = Real(1e-10),
int lim = 24) {
assert(eps > Real{});
assert(lim >= 0);
if (l == r) {
return Real{};
}
bool rev = r < l;
if (rev) std::swap(l, r);
// Keep the usual arithmetic when finite; split before adding on overflow.
auto midpoint = [](Real a, Real b) {
Real s = a + b;
return std::isfinite(s) ? s / Real(2) : a / Real(2) + b / Real(2);
};
auto sim = [](Real a, Real b, Real fa, Real fm, Real fb) {
Real wid = b - a;
Real val = wid * (fa + Real(4) * fm + fb) / Real(6);
if (std::isfinite(val)) return val;
Real avg = std::lerp(std::midpoint(fa, fb), fm, Real(2) / Real(3));
return std::isfinite(wid) ? wid * avg
: ((b / Real(2) - a / Real(2)) * avg) * Real(2);
};
Real mdd = midpoint(l, r);
Real fl = std::invoke(fun, l);
Real fm = std::invoke(fun, mdd);
Real fr = std::invoke(fun, r);
Real sum = sim(l, r, fl, fm, fr);
auto sol = [&](auto &dfs, Real a, Real b, Real mid, Real fa, Real fc,
Real fb, Real est, Real tol, int dep) -> Real {
if (mid == a || mid == b) return (b - a) * std::midpoint(fa, fb);
Real lm = midpoint(a, mid);
Real rm = midpoint(mid, b);
Real flm = std::invoke(fun, lm);
Real frm = std::invoke(fun, rm);
Real sl = sim(a, mid, fa, flm, fc);
Real sr = sim(mid, b, fc, frm, fb);
Real dif = sl + sr - est;
if (dep == 0 || std::abs(dif) <= Real(15) * tol) {
return sl + sr + dif / Real(15);
}
return dfs(dfs, a, mid, lm, fa, flm, fc, sl, tol / Real(2), dep - 1) +
dfs(dfs, mid, b, rm, fc, frm, fb, sr, tol / Real(2), dep - 1);
};
Real res = sol(sol, l, r, mdd, fl, fm, fr, sum, eps, lim);
return rev ? -res : res;
}
} // namespace noya
#endif // NOYA_ADAPTIVE_SIMPSON_HPP
/// @complexity Time: O(2^d) function evaluations in the worst case, where d is lim.
/// Space: O(d) recursion stack.
#include <cassert>
#include <cmath>
#include <functional>
#include <numeric>
#include <utility>
namespace noya {
/// @brief Numerically integrate a smooth real function on [l, r] with
/// adaptive Simpson subdivision. eps targets the estimated absolute error;
/// recursion also stops at lim or floating-point resolution, without a certified bound.
template <class Function, class Real = long double>
Real adaptive_simpson(Function fun, Real l, Real r, Real eps = Real(1e-10),
int lim = 24) {
assert(eps > Real{});
assert(lim >= 0);
if (l == r) {
return Real{};
}
bool rev = r < l;
if (rev) std::swap(l, r);
// Keep the usual arithmetic when finite; split before adding on overflow.
auto midpoint = [](Real a, Real b) {
Real s = a + b;
return std::isfinite(s) ? s / Real(2) : a / Real(2) + b / Real(2);
};
auto sim = [](Real a, Real b, Real fa, Real fm, Real fb) {
Real wid = b - a;
Real val = wid * (fa + Real(4) * fm + fb) / Real(6);
if (std::isfinite(val)) return val;
Real avg = std::lerp(std::midpoint(fa, fb), fm, Real(2) / Real(3));
return std::isfinite(wid) ? wid * avg
: ((b / Real(2) - a / Real(2)) * avg) * Real(2);
};
Real mdd = midpoint(l, r);
Real fl = std::invoke(fun, l);
Real fm = std::invoke(fun, mdd);
Real fr = std::invoke(fun, r);
Real sum = sim(l, r, fl, fm, fr);
auto sol = [&](auto &dfs, Real a, Real b, Real mid, Real fa, Real fc,
Real fb, Real est, Real tol, int dep) -> Real {
if (mid == a || mid == b) return (b - a) * std::midpoint(fa, fb);
Real lm = midpoint(a, mid);
Real rm = midpoint(mid, b);
Real flm = std::invoke(fun, lm);
Real frm = std::invoke(fun, rm);
Real sl = sim(a, mid, fa, flm, fc);
Real sr = sim(mid, b, fc, frm, fb);
Real dif = sl + sr - est;
if (dep == 0 || std::abs(dif) <= Real(15) * tol) {
return sl + sr + dif / Real(15);
}
return dfs(dfs, a, mid, lm, fa, flm, fc, sl, tol / Real(2), dep - 1) +
dfs(dfs, mid, b, rm, fc, frm, fb, sr, tol / Real(2), dep - 1);
};
Real res = sol(sol, l, r, mdd, fl, fm, fr, sum, eps, lim);
return rev ? -res : res;
}
} // namespace noya
// NOYA_ADAPTIVE_SIMPSON_HPP