From 00245b1d61cc6f9faac7834b94138f1dfec4172c Mon Sep 17 00:00:00 2001
From: i1011 <56438495+i1011@users.noreply.github.com>
Date: Sun, 15 Oct 2023 21:39:37 +0800
Subject: [PATCH] Update AdaptiveSimpson.cpp

---
 codebook/9_Else/AdaptiveSimpson.cpp | 58 +++++++++++++++--------------
 1 file changed, 30 insertions(+), 28 deletions(-)

diff --git a/codebook/9_Else/AdaptiveSimpson.cpp b/codebook/9_Else/AdaptiveSimpson.cpp
index 54328d20..49a8357a 100644
--- a/codebook/9_Else/AdaptiveSimpson.cpp
+++ b/codebook/9_Else/AdaptiveSimpson.cpp
@@ -1,28 +1,30 @@
-using F_t = function<double(double)>;
-pdd simpson(const F_t &f, double l, double r,
-  double fl, double fr, double fm = nan("")) {
-  if (isnan(fm)) fm = f((l + r) / 2);
-  return {fm, (r - l) / 6 * (fl + 4 * fm + fr)};
-}
-double simpson_ada(const F_t &f, double l, double r,
-  double fl, double fm, double fr, double eps) {
-  double m = (l + r) / 2,
-         s = simpson(f, l, r, fl, fr, fm).second;
-  auto [flm, sl] = simpson(f, l, m, fl, fm);
-  auto [fmr, sr] = simpson(f, m, r, fm, fr);
-  double delta = sl + sr - s;
-  if (abs(delta) <= 15 * eps)
-    return sl + sr + delta / 15;
-  return simpson_ada(f, l, m, fl, flm, fm, eps / 2) +
-    simpson_ada(f, m, r, fm, fmr, fr, eps / 2);
-}
-double simpson_ada(const F_t &f, double l, double r) {
-  return simpson_ada(
-    f, l, r, f(l), f((l + r) / 2), f(r), 1e-9 / 7122);
-}
-double simpson_ada2(const F_t &f, double l, double r) {
-    double h = (r - l) / 7122, s = 0;
-  for (int i = 0; i < 7122; ++i, l += h)
-    s += simpson_ada(f, l, l + h);
-  return s;
-}
+template<typename Func, typename d = double>
+struct Simpson {
+  using pdd = pair<d, d>;
+  Func f;
+  pdd mix(pdd l, pdd r, optional<d> fm = {}) {
+    d h = (r.X - l.X) / 2, v = fm.value_or(f(l.X + h));
+    return {v, h / 3 * (l.Y + 4 * v + r.Y)};
+  }
+  d eval(pdd l, pdd r, d fm, d eps) {
+    pdd m((l.X + r.X) / 2, fm);
+    d s = mix(l, r, fm).second;
+    auto [flm, sl] = mix(l, m);
+    auto [fmr, sr] = mix(m, r);
+    d delta = sl + sr - s;
+    if (abs(delta) <= 15 * eps) return sl + sr + delta / 15;
+    return eval(l, m, flm, eps / 2) +
+      eval(m, r, fmr, eps / 2);
+  }
+  d eval(d l, d r, d eps) {
+    return eval({l, f(l)}, {r, f(r)}, f((l + r) / 2), eps);
+  }
+  d eval2(d l, d r, d eps, int k = 997) {
+    d h = (r - l) / k, s = 0;
+    for (int i = 0; i < k; ++i, l += h)
+      s += eval(l, l + h, eps / k);
+    return s;
+  }
+};
+template<typename Func>
+Simpson<Func> make_simpson(Func f) { return {f}; }