自适应辛普森（算法简要 + 模板）

简述：

对于普通的二次函数有f(x)=ax2+bx+cf(x) = ax ^ 2 + bx + cf(x)=ax2+bx+c，原函数F(x)=ax33+bx22+cxF(x) = \frac{a x^3}{3} + \frac{bx ^2}{2} + cxF(x)=3ax3+2bx2+cx。

有
⎰lrf(x)=F(r)−F(l)=a3(r−l)3+b2(r−l)2+c(r−l)=(r−l)6(2a(l2+r2+lr)+3b(l+r)+c)=(r−l)6((al2+bl+c)+(ar2+br+c)+4(a((l+r)2)2+b(l+r)2+c))=(r−l)6(f(l)+f(r)+4f((l+r)2))\lmoustache_{l} ^{r} f(x) = F(r) - F(l) \\ = \frac{a}{3}(r - l) ^3 + \frac{b}{2} (r - l) ^2 + c(r - l)\\ = \frac{(r - l)}{6} (2a(l ^ 2 + r ^2 + lr) + 3b(l + r) + c)\\ = \frac{(r - l)}{6} ((al ^2 + bl + c) +( a r ^ 2 + br + c) + 4(a(\frac{(l + r)}{2}) ^ 2 + b\frac{(l + r)}{2} + c))\\ = \frac{(r - l)}{6} (f(l) + f(r) + 4f(\frac{(l + r)}{2}))\\ ⎰lrf(x)=F(r)−F(l)=3a(r−l)3+2b(r−l)2+c(r−l)=6(r−l)(2a(l2+r2+lr)+3b(l+r)+c)=6(r−l)((al2+bl+c)+(ar2+br+c)+4(a(2(l+r))2+b2(l+r)+c))=6(r−l)(f(l)+f(r)+4f(2(l+r)))
这是普通的辛普森推导，自适应辛普森就是在普通的辛普森上加了一个评估步骤，我们通过把一个区间二分，如果这个区间整体求得的值与一半加一半单独求得的值是较为接近的这个时候我们就可以认为我们已经求得正确答案。

P4525 【模板】自适应辛普森法1

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-6;double a, b, c, d, l, r;double f(double x) {return (c * x + d) / (a * x + b);
}double sim(double l, double r) {return (r - l) * (f(l) + f(r) + 4.0 * f((l + r) / 2.0)) / 6.0;
}double asr(double l, double r, double eps, double ans) {double mid = (l + r) / 2;double ansl = sim(l, mid), ansr = sim(mid, r);if(fabs(ansl + ansr - ans) < 15 * eps) return ansl + ansr + (ansl + ansr - ans) / 15.0;return asr(l, mid, eps / 2.0, ansl) + asr(mid, r, eps / 2.0, ansr);
}int main() {// freopen("in.txt", "r", stdin);   // freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &l, &r);printf("%.6f\n", asr(l, r, eps, sim(l, r)));return 0;
}

Ellipse

显然有y=b×1−(xa)2y = b \times \sqrt{1 - (\frac{x}{a}) ^2}y=b×1−(ax)2，所以只要简单的带入求解即可。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int N = 1e3 + 10, mod = 6;double a, b, l, r;double f(double x) {return 2 * b * (sqrt(1 - (x * x) / (a * a)));
}double sim(double l, double r) {return (r - l) * (f(l) + f(r) + 4.0 * f((l + r) / 2.0)) / 6.0;
}double asr(double l, double r, double eps, double ans) {double mid = (l + r) / 2;double ansl = sim(l, mid), ansr = sim(mid, r);if(fabs(ansl + ansr - ans) < eps) return ansl + ansr;return asr(l, mid, eps / 2.0, ansl) + asr(mid, r, eps / 2.0, ansr);
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int T; scanf("%d", &T);while(T--) {scanf("%lf %lf %lf %lf", &a, &b, &l, &r);printf("%.3f\n", asr(l, r, 1e-5, sim(l, r)));}return 0;
}

The area

两个函数积分求差，也就是分别算出两个积分，然后相减即可。

对于二次函数的式子，这里直接用了拉格朗日插值来求了，没去推导公式。

对于一次函数的式子，只要简单求解一下即可。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-6;double x[5], y[5];double f1(double X) {double ans = 0;for(int i = 1; i <= 3; i++) {double a = 1, b = 1;for(int j = 1; j <= 3; j++) {if(i == j) continue;a *= X - x[j];b *= x[i] - x[j];}ans += y[i] * a / b;}return ans;
}double sim1(double l, double r) {return (r - l) * (f1(l) + f1(r) + f1((l + r) / 2.0) * 4.0) / 6.0;
}double asr1(double l, double r, double eps, double ans) {double mid = (l + r) / 2.0;double ansl = sim1(l, mid), ansr = sim1(mid, r);if(fabs(ansl + ansr - ans) < 15.0 * eps) return ansl + ansr + (ansl + ansr - ans) / 25.0;return asr1(l, mid, eps / 2.0, ansl) + asr1(mid, r, eps / 2.0, ansr);
}double f2(double X) {return (y[3] - y[2]) / (x[3] - x[2]) * (X - x[3]) + y[3];
}double sim2(double l, double r) {return (r - l) * (f2(l) + f2(r) + f2((l + r) / 2.0) * 4.0) / 6.0;
}double asr2(double l, double r, double eps, double ans) {double mid = (l + r) / 2;double ansl = sim2(l, mid), ansr = sim2(mid, r);if(fabs(ansl + ansr - ans) < 15.0 * eps) return ansl + ansr + (ansl + ansr - ans) / 15.0;return asr2(l, mid, eps / 2.0, ansl) + asr2(mid, r, eps / 2.0, ansr);
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int T; scanf("%d", &T);while(T--) {for(int i = 1; i <= 3; i++) {scanf("%lf %lf", &x[i], &y[i]);}double ans = asr1(x[2], x[3], eps, sim1(x[2], x[3])) - asr2(x[2], x[3], eps, sim2(x[2], x[3]));printf("%.2f\n", ans);}return 0;
}