adaptive_simpson.hpp¶

SECTIONOptimization INCLUDEnoya/adaptive_simpson.hpp

按误差估计自适应细分区间，数值计算平滑函数的定积分。

Complexity: Time: O(2^d) function evaluations in the worst case, where d is lim. Space: O(d) recursion stack.

Implementation¶

ReferenceSourceExpanded

当前头文件，省略 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>

namespace noya {

/// @brief Numerically integrate a smooth real function on [l, r] with
/// adaptive Simpson subdivision and an absolute error target.
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{};
  }
  if (r < l) {
    return -adaptive_simpson<Function, Real>(fun, r, l, eps, lim);
  }

  auto sim = [](Real a, Real b, Real fa, Real fm, Real fb) {
    return (b - a) * (fa + Real(4) * fm + fb) / Real(6);
  };
  Real mdd = (l + r) / Real(2);
  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 fa, Real fc, Real fb, Real est,
                 Real tol, int dep) -> Real {
    Real mid = (a + b) / Real(2);
    Real lm = (a + mid) / Real(2);
    Real rm = (mid + b) / Real(2);
    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, fa, flm, fc, sl, tol / Real(2), dep - 1) +
           dfs(dfs, mid, b, fc, frm, fb, sr, tol / Real(2), dep - 1);
  };
  return sol(sol, l, r, fl, fm, fr, sum, eps, lim);
}

} // namespace noya

#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>

namespace noya {

/// @brief Numerically integrate a smooth real function on [l, r] with
/// adaptive Simpson subdivision and an absolute error target.
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{};
  }
  if (r < l) {
    return -adaptive_simpson<Function, Real>(fun, r, l, eps, lim);
  }

  auto sim = [](Real a, Real b, Real fa, Real fm, Real fb) {
    return (b - a) * (fa + Real(4) * fm + fb) / Real(6);
  };
  Real mdd = (l + r) / Real(2);
  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 fa, Real fc, Real fb, Real est,
                 Real tol, int dep) -> Real {
    Real mid = (a + b) / Real(2);
    Real lm = (a + mid) / Real(2);
    Real rm = (mid + b) / Real(2);
    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, fa, flm, fc, sl, tol / Real(2), dep - 1) +
           dfs(dfs, mid, b, fc, frm, fb, sr, tol / Real(2), dep - 1);
  };
  return sol(sol, l, r, fl, fm, fr, sum, eps, lim);
}

} // namespace noya

#endif // NOYA_ADAPTIVE_SIMPSON_HPP

#include <cassert>
#include <cmath>
#include <functional>

/// @complexity Time: O(2^d) function evaluations in the worst case, where d is lim.
/// Space: O(d) recursion stack.

namespace noya {

/// @brief Numerically integrate a smooth real function on [l, r] with
/// adaptive Simpson subdivision and an absolute error target.
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{};
  }
  if (r < l) {
    return -adaptive_simpson<Function, Real>(fun, r, l, eps, lim);
  }

  auto sim = [](Real a, Real b, Real fa, Real fm, Real fb) {
    return (b - a) * (fa + Real(4) * fm + fb) / Real(6);
  };
  Real mdd = (l + r) / Real(2);
  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 fa, Real fc, Real fb, Real est,
                 Real tol, int dep) -> Real {
    Real mid = (a + b) / Real(2);
    Real lm = (a + mid) / Real(2);
    Real rm = (mid + b) / Real(2);
    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, fa, flm, fc, sl, tol / Real(2), dep - 1) +
           dfs(dfs, mid, b, fc, frm, fb, sr, tol / Real(2), dep - 1);
  };
  return sol(sol, l, r, fl, fm, fr, sum, eps, lim);
}

} // namespace noya