自适应辛普森法
题意：
试计算积分
\[ \displaystyle{\int_L^R\frac{cx+d}{ax+b}\,\mathrm{d}x} \]
结果保留至小数点后 \(6\) 位。
数据保证计算过程中分母不为 \(0\) 且积分能够收敛。
输入格式：
一行，包含 \(6\) 个实数 \(a,b,c,d,L,R\)。
输出格式：
一行，积分值，保留至小数点后 \(6\) 位。
数据范围：
\(a,b,c,d\in[-10,10]\)，\(-100\le L<R\le 100\) 且 \(R-L\ge1\)。
先给出 Simpson 公式 \[ \int_a^b f(x) \,\mathrm{d} x \,\approx\,(b-a)\cdot \frac{f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)}{6} \] 下面给出推导和介绍
Simpson 公式的原理在于用二次函数去拟合区间，并计算定积分。
显然，区间划分的越多，精确度也就相对越高，不过相应的时间复杂度也会增大。
不过我们先来看看 Simpson 公式是怎么推导的。 \[ \begin{aligned} \int_a^b f(x) \,\mathrm{d} x &\approx \int_a^b (A x^2+B x+C) \,\mathrm{d} x \\[6pt]&=\frac{A}{3}\left(b^3-a^3\right)+\frac{B}{2}\left(b^2-a^2\right)+C(a-b) \\[6pt]&=\frac{(b-a)}{6}\left(2 A\left(b^2+a b+a^2\right)+3 B(b+a)+6 C\right) \\[6pt]&=\frac{(b-a)}{6}\left(2 A b^2+2 A a b+2 A a^2+3 B b+3 B a+6 C\right) \\[6pt]&=\frac{(b-a)}{6}\left(A a^2+B a+C+A b^2+B b+C+4 A\left(\frac{a+b}{2}\right)^2+4 B\left(\frac{a+b}{2}\right)+4 C\right) \\[6pt]&=\frac{(b-a)}{6}\left(f(a)+f(b)+4 f\left(\frac{a+b}{2}\right)\right) \end{aligned} \] 这里第一行用到了积分的小知识，可以看我之前写的这篇不定积分基本公式
接着我们来讲本题用到的自适应辛普森法。
刚刚我们提到了 Simpson 公式需要根据精度进行分段拟合，而自适应辛普森法就是实现这个操作的算法
算法原理也很简单，就是把原来的区间二分递归，如果精度不够就继续递归，否则直接退出。
时间复杂度取决于原函数 \(f(x)\) 的复杂程度。
代码：
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)())
const double eps = 1e-12;
double a,b,c,d,l,r;
double f(double x) { return (c * x + d) / (a * x + b); }
double simpson(double l, double r)
{
double mid = (l + r) / 2;
return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double calc(double l, double r,double ans)
{
double mid = (l + r) / 2;
double L = simpson(l, mid), R = simpson(mid, r);
if(fabs(L + R - ans) <= eps) return L + R + (L + R - ans);
return calc(l, mid, L) + calc(mid, r, R);
}
signed main()
{
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
// freopen("check.in","r",stdin);
// freopen("check.out","w",stdout);
cin >> a >> b >> c >> d >> l >> r;
cout << fixed << setprecision(6) << calc(l, r, simpson(l, r)) << '\n';
return 0;
}参考文献：
