# 自适应辛普森法

Author: Fanyi Pu

Published: 2019-02-04

Canonical: <https://pufanyi.com/blog/oi-icpc/notes/simpson>

辛普森公式的推导与自适应积分实现。

并没想打算开篇小小结，只是随便写写（感觉开小小结的成本太高了啊……）

感觉还是挺妙的。以前听[zyy 大佬](https://blog.csdn.net/zhouyuyang233)讲课的时候说过，函数的积分用二次函数拟合一下。

加入我们要算 $\int_l^r F(x)\mathrm{d}x$，我们用二次函数 $f(x)=ax^2+bx+c$ 来拟合。

我们发现：

$$
\begin{aligned}
\int_l^r f(x)\mathrm{d}x&=\frac{a}{3}r^3+\frac{b}{2}r^2+cr-\frac{a}{3}l^3-\frac{b}{2}l^2-cl\\
&=\frac{a}{3}\left(r-l\right)\left(r^2+l^2+lr\right)+\frac{b}{2}\left(r-l\right)\left(r+l\right)+c\left(r-l\right)\\
&=\frac{r-l}{6}\left(2ar^2+2al^2+2alr+3br+3bl+6c\right)\\
&=\frac{r-l}{6}\left(ar^2+br+c+al^2+bl+c+ar^2+al^2+2alr+2br+2bl+4c\right)\\
&=\frac{r-l}{6}\left(f(l)+f(r)+a(l+r)^2+2b(l+r)+4c\right)\\
&=\frac{r-l}{6}\left(f(l)+f(r)+4\left(a\left(\frac{l+r}{2}\right)^2+b\left(\frac{l+r}{2}\right)+c\right)\right)\\
&=\frac{r-l}{6}\left(f(l)+f(r)+4f\left(\frac{l+r}{2}\right)\right)
\end{aligned}
$$

于是我们有：

$$
\int_l^r F(x)\mathrm{d}x\approx\frac{r-l}{6}\left(F(l)+F(r)+4F\left(\frac{l+r}{2}\right)\right)
$$

然后就可以拟合辣！

```cpp
inline double simpson(double l, double r)
{
	return (F(l) + F(r) + 4 * F((l + r) / 2)) * (r - l) / 6;
}
```

然而就是要么 WA（炸精度）要么 T……

有一种叫自适应辛普森的方法解决这个问题，我们用辛普森计算 $\int_l^rF(x)\mathrm{d}x$，$\int_l^{(l+r)/2}F(x)\mathrm{d}x$，$\int_{(l+r)/2}^rF(x)\mathrm{d}x$ 的值，如果

$$
\left|\int_l^rF(x)\mathrm{d}x-\left(\int_l^{(l+r)/2}F(x)\mathrm{d}x+\int_{(l+r)/2}^rF(x)\mathrm{d}x\right)\right|\le \mathrm{eps}
$$

那么直接返回值，否则递归下去算。

upd：忽然发现我学了个假的 `simpson`。其实 `eps` 应该是相对误差，所以每次递归下去时 $\mathrm{eps}$ 都应该除 $2$。听说最后判定时应该是 $\le 15\cdot \mathrm{eps}$，虽然我并不知道是为什么。话说这和 eps 乘上 15 有什么区别吗

这样复杂度就是 $\mathcal{O}\left(\text{玄学}\right)$，交之前洗把脸就行了。

```cpp
inline double asr(double l, double r, double ans, double eps)
{
	double mid = (l + r) * .5;
	double ll = simpson(l, mid), rr = simpson(mid, r);
	if(fabs(ll + rr - ans) < 15 * eps)
		return ll + rr;
	return asr(l, mid, ll, eps * .5) + asr(mid, r, rr, eps * .5);
}
```