\(\rm{0x01}\) 前言
首先闡明一點，自適應辛普森演算法(\(\rm{Adaptive ~Simpson’s~ rule}\) )是一類近似演算法(\(\rm{Approximation ~algorithm}\))，主要用於求較難反導的函式的積分。大概在資訊計算的時候中很常用？
其思想是利用二次函式來不斷擬合(\(\rm{Overfitting}\))所求曲線，而所謂的\(Adapative\)（自適應）則是用於優化時間複雜度的方法。
嗝…總之…比較簡單？
表示看了兩篇外國學者的論文，感覺好像思路已經比較清晰的樣子。
\(\rm{0x02}\) \(\mathcal{Simpson~Formula}\) 辛普森公式
稍等，這個跟演算法的關係不大，主要是公式：\[\rm{\int _{l}^{r} f(x) \rm{dx}\approx \frac{f(r) +f(l)+ 4 \cdot f(m)}{6}} \,\]
事實上吧，求積分的話，大多數都是直接套辛普森公式的。並且這個公式是廣泛適用的……只不過誤差有點兒人大\(233\)
這其實是我們用二次函式不斷擬合的結果，推導過程大致如下\(^{[1]}\)：
因為 \(g(x)\) 是用來擬合 \(f(x)\) 的，所以有：
\(\int\limits_L^Rf(x)dx\approx \int\limits_L^Rax^2+bx+c\space dx\)
求 \(g(x)\) 的不定積分為：
\(\int g(x)dx=\frac13ax^3+\frac12bx^2+cx+C\)
然後再帶入 RR 和 LL ：
\(\int\limits_L^Rf(x)dx=\frac13a(R^3-L^3)+\frac12b(R^2-L^2)+c(R-L)\)
然後提公因式，原式為：
\(\frac{R-L}{6}(2a(R^2+LR+L^2)+3b(R+L)+6c)6R−L(2a(R2+LR+L2)+3b(R+L)+6c)\)
把裡面展開來：
\(\frac{R-L}{6}(2aR^2+2aLR+2aL^2+3bR+3bL+6c)\)
重新整理一下式子：
\(\frac{R-L}{6}((aR^2+bR+c)+(aL^2+bL+c)+aR^2+aL^2+2aLR+2bR+2bL+4c)\)
再整理：
\(\frac{R-L}{6}((aR^2+bR+c)+(aL^2+bL+c)+4(a(\frac{R+L}{2})^2+b(\frac{R+L}{2})+c))\)
代換可得：
\(aR^2+bR+c\approx f(R)\)
\(aL^2+bL+c\approx f(L)\)
\(4(a(\frac{R+L}{2})^2+b(\frac{R+L}{2})+c)\approx 4f(\frac{R+L}{2})\)
把這三個式子帶回去， 最後我們就得到了\(\rm{Simpson}\) 積分公式：\(\int\limits_L^Rf(x)dx \approx \frac{R-L}{6}(f(L)+4f(\frac{L+R}{2})+f(R)))\)
於是我們就得到了所謂的\(\rm{Simpson~Fomula}\)。但事實上，對於一段“跌宕起伏”的函式，我們還是無法準確地去用一個二次函式的擬合來刻畫。於是——
\(\rm{0x03}\) \(\mathcal{Adaptive~Simpson~Rule}\) 自適應辛普森法
我們考慮，如果把區間們稍微劃分地更細一點，那是不是會好呢？答案是一定的。那麼我們可以考慮定向二分。但是……定向二分始終存在一個問題，就是它太笨了，笨到如果\([l_1,r_1]\)已經滿足精度要求了，我們卻還要一直分；笨到\([l_2,r_2]\)分的次數根本不夠用——但我們並不可以得到反饋。
於是考慮自適應。
所謂自適應，說的直白點，無非就是需要多分治幾次的地方，多分治幾次；不需要的則可以少分治幾次
你會發現，其實他節約的就是一個點——時間效率。
舉個例子\(^{[2]}\)：
比如有這麼個函式\(\rm{f (x) = 13(x − x^2)e^{\frac{−3x}{2}} }\)，我們要求\(\int_{0}^{4}f(x) \rm{~dx}\) 並要求精度誤差在\(1e-5\) 以內。而我們有兩種方法去解決：
- 以固定的比例以及約束二分。
- 運用自適應策略分配
那麼我們首先要知道他真正的\(value:\)
看起來好像海星？然後我們用兩種方法都試一試：
首先是自適應法，我們發現最後只需要求\(20\)段區間。表中的\(a_k \& b_k\)表示左右區間，\(S(l,r)\)表示\([l,r]\)內、運用\(0x01\)中的公式計算的，\(\rm{Error~Bound}\)表示誤差界，\(\rm{Tolerance}\)表示計算時需要的誤差(程式設計時會講)。
那麼最後算出來的值是$ −1.54878823413$ ，與真實值誤差為$0.00000013840 $，一共呼叫了\(79\)次函式估值（留個坑，後文會講）。
那麼繪製出來的函式影像大概長這樣：
好像很流暢？\(233\)
那麼第二種方法是定值分段。我們考慮分成區間\([0,4]\)分為長度為\(0.03125\)的\(128\)段，並運用\(0x01\)的\(Formula\)，最後得出的結果為\(−1.54878844029\)，誤差為\(0.00000006776\)……
好像是第一個誤差的二分之一？聽起來好像誤差小了，但是卻需要\(257\)次函式估值的呼叫……相比之下，我們可以獲得更優越的效能，而那誤差也是不需要考慮的啦的啦～
但是比起\(1e-5\)精度來說……這波穩賺啊\(233\)
\(\rm{0x04}\) \(\mathcal{About~Code~}\) 程式碼實現
首先是\(LuoguP4525\)的暴力解法：
#include <cmath>
#include <cstdio>
#include <iostream>
#define mid (l + r) / 2
using namespace std ;
const double Eps = 5e-8 ;
double Ans, A, B, C, D, L, R ;
inline double F(double x){
return (C * x + D) / (A * x + B) ;
}
inline double Simp_calc(double l, double r){
return (r - l) / 6 * (F(l) + F(r) + 4 * F(mid)) ;
}
inline double do_divide(double l, double r){
double Lv, Rv, v, t ;
Lv = Simp_calc(l, mid), Rv = Simp_calc(mid, r), v = Simp_calc(l, r) ;
if (fabs(t = (Lv + Rv - v)) <= Eps) return Lv + Rv ; return do_divide(l, mid) + do_divide(mid, r) ;
}
int main(){
scanf("%lf%lf%lf%lf%lf%lf", &A, &B, &C, &D, &L, &R) ;
Ans = do_divide(L, R) ; printf("%.6lf", Ans) ; return 0 ;
}這…嚴格意義上講，不能算是自適應辛普森法——更準確地說，四不像，或者“東施效顰”之類的，都挺合適。這是我在初學這一塊兒內容時的寫法，他不嚴格正確，但是……他對的很？
至於進化版\(LuoguP4526\)，也是完全可以violently艹過去的：
#include <cmath>
#include <cstdio>
#include <iostream>
#define mid (l + r) / 2
using namespace std ;
const double Eps = 5e-8 ;
double Ans, A, B, C, D, L, R ;
inline double F(double x){
return pow(x, A / x - x) ;
}
inline double Simp_calc(double l, double r){
return (r - l) / 6 * (F(l) + F(r) + 4 * F(mid)) ;
}
inline double do_divide(double l, double r, double Lans){
double Lv, Rv, v, t ; Lv = Simp_calc(l, mid), Rv = Simp_calc(mid, r), v = Lans ;
if (fabs(t = (Lv + Rv - v)) <= Eps) return Lv + Rv ; return do_divide(l, mid, Lv) + do_divide(mid, r, Rv) ;
}
int main(){
scanf("%lf", &A) ; L = Eps, R = 23.3 ;
if (A < 0) { puts("orz"), cout << endl ; return 0 ;}
Ans = do_divide(L, R, Simp_calc(L, R)) ; printf("%.5lf", Ans) ; return 0 ;
}但是，其真正的實現應該是這樣：
#include <cmath>
#include <cstdio>
#include <iostream>
#define v Lans
#define hps eps * 0.5
#define mid (l + r) * 0.5
using namespace std ;
double Ans, A, L, R ;
const double Eps = 1e-7 ;
const double Liu = 1.0 / 6 ;
inline double F(double x){
return pow(x, A / x - x) ;
}
inline double Simp_calc(double l, double r){
return (r - l) * Liu * (F(l) + F(r) + 4.0 * F(mid)) ;
}
inline double do_divide(double l, double r, double Lans, double eps){
double Lv, Rv, t ;
Lv = Simp_calc(l, mid), Rv = Simp_calc(mid, r) ;
if (fabs(t = (Lv + Rv - v)) <= eps * 15) return Lv + Rv + t / 15 ;
return do_divide(l, mid, Lv, hps) + do_divide(mid, r, Rv, hps) ; //據說eps×15來自於Wiki……
}
int main(){
scanf("%lf", &A) ; L = Eps, R = 30 ;
if (A < 0) { puts("orz"), cout << endl ; return 0 ;}
Ans = do_divide(L, R, Simp_calc(L, R), Eps) ; printf("%.5lf", Ans) ; return 0 ;
}然後我們發現這就可以跟文章上面\(0x03\)中的例子呼應了：每次分治時計算兩次，總共分治了\(39\)次，最終一共計算了\(78+1=79\)次，而二分則是一棵有\(128\)個葉子節點的遞迴樹，總共計算了\(256 +1=257\)次。
好的，終於要扯正題了。演算法的實現其實簡單，我們用擬合性演算法不斷\(check\&calc\), 而\(check\)的方式也很簡單，只需要判斷一下兩段子區間的函式值之和與整個區間的函式值之和的差值是否在精度要求範圍之內，之後如果滿足精度誤差就直接\(return\)，否則對於這段區間繼續遞迴下去。
而這個地方有個要求，就是對於\(eps\)，你需要不斷\(half\)他，原因很簡單，對於一整段區間\(U\)，要求他的返回值的\(|eps(U)| \leq k\)的話，那麼對於其作為子集兩個連續區間\(A,B\)，當\(A \bigcup B = U\)時，必須要滿足\(|eps(A)| \leq \frac{k}{2}, |eps(B)| \leq \frac{k}{2}\)，才能保證\(|eps(U) = eps(A) + eps(B)| \leq k\)，所以要：
#define hps eps * 0.5
....................
return do_divide(l, mid, Lv, hps) + do_divide(mid, r, Rv, hps) ; 好了，唯一的問題在於有句話迷的很：
if (fabs(t = (Lv + Rv - v)) <= eps * 15) return Lv + Rv + t / 15 ; 
這個\(\leq 15 \cdot eps\)是什麼意思？
好吧，筆者也不清楚，但是有個結論是長這樣的：
什麼？你說推什麼倒？小小年紀整天想著推倒學姐可不好啊\(233\)
什麼？你還想看推倒的過程？嘖嘖嘖，左轉知乎感情專區蟹蟹～
好的，以上兩段是瞎扯的，推導過程十分詭異，大家自己食用好了……
於是結束，撒花花…（不是撒博主x）
\(\rm{Referance}\)
\([1]\) :\(NaCly\_Fish\)的\(blog\) \(^{^{[\nearrow]}}\)
\([2]\) :某篇論文 \(^{^{[\nearrow]}}\)
對於第二篇\(refer\)，借鑑的比較多（讀書人的事……咳咳……），但是改了一個資料，就是\(81 \to 79\)，因為程式碼需要\(233\)
\(Ps:\)史上最不正經的\(reference\)誕生了……