前言

蒟蒻的计算几何简直太垃圾了，什么都不会，所以要从最简单的东西开始......要是出现错误，请在下方评论区指出。

为了让某一天 或者每天 犯傻的 $CXY07$ 可以看得懂，一切都用尽我所能用最 sb 显然的语言解释，巨佬们不要嫌我啰嗦 $qwq$

正题

参考博客

如果我们需要求解这样一个东西：

$$\int_l^r f(x)\,dx$$

我们显然有两种方法：

1、直接极限搞，把这段区间分割为无数个小区间，就有：

$$\int_l^r f(x)\,dx = \lim_{n \to \infty} \sum_{i=1}^n \dfrac{r-l}{n} \times f(l+\dfrac{r-l}{n}\times i)$$
2、牛顿-莱布尼茨公式法，令 $F'(x)=f(x)$，那么显然有：

$$\int_l^r f(x)\,dx = F(r)-F(l)$$

这两种方法都有一定的适应范围，但对于一些两种方法都不适用的呢，我们就需要自适应辛普森法......

辛普森公式

辛普森公式其实思路很简单，对于一些我们不是很好积分的函数，直接拿跟这个函数差不多的二次函数近似它，然后去积二次函数

假设 $f(x)$ 用二次函数近似后是 $g(x)=ax^2+bx+c$

那我们有这样一个东西：

$\int_l^r f(x)\,dx \approx \int_l^r g(x)\,dx$

$=\int_l^r (ax^2+bx+c)\,dx$

$=\dfrac{a}{3} \times (r^3-l^3) + \dfrac{b}{2} \times (r^2-l^2) + c \times (r - l)$

$=\dfrac{2a\times (r^3-l^3) + 3b \times (r^2-l^2) + 6c \times (r-l)}{6}$

$=\dfrac{2a\times (r-l)(r^2+lr+l^2) + 3b \times (r-l)(r+l) + 6c \times (r-l)}{6}$

$=\dfrac{(r-l)(2a\times (r^2+lr+l^2) + 3b \times (r+l) + 6c)}{6}$

$=\dfrac{(r-l)}{6} \times (2ar^2+2alr+2al^2+3br+3bl+6c)$

$=\dfrac{(r-l)}{6} \times (al^2+bl+c+ar^2+br+c+ar^2+2alr+al^2+2br+2bl+4c)$

$=\dfrac{(r-l)}{6} \times (g(l)+g(r)+ar^2+2alr+al^2+2br+2bl+4c)$

$=\dfrac{(r-l)}{6} \times (g(l)+g(r)+a(l+r)^2+2b(l+r)+4c)$

$=\dfrac{(r-l)}{6} \times (g(l)+g(r)+4(a(\dfrac{l+r}{2})^2+b(\dfrac{l+r}{2})+c))$

$=\dfrac{(r-l)}{6} \times (g(l)+g(r)+4g(\dfrac{l+r}{2}))$

$\approx \dfrac{(r-l)(f(l)+f(r)+4f(\dfrac{l+r}{2}))}{6}$

我们就得出了这个公式：

$$\int_l^r f(x)\,dx \approx \dfrac{(r-l)(f(l)+f(r)+4f(\dfrac{l+r}{2}))}{6}$$

代码大概长这样：

inline double simpson(double l,double r) {
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}

显然这东西误差很大，如果是一个长得不那么二次函数的函数，那就得凉凉，但是可以发现 $(r-l)$ 越小，这东西应该会越接近 感性理解

那么我们现在就需要处理精度问题了

自适应辛普森法（ASR）

大概思路就是把 $[l,r]$ 多去拆分为几个小区间，每个小区间的 $[l,r]$ 显然更小，那么就应该更接近答案，这个东西递归就行

但是分少了精度不够，分多了 $TLE$，我们就需要一个方法来让程序自己判断是否需要继续递归

我们设 $[l,r]$ 区间内直接用辛普森公式的答案为 $ans$，令 $mid=\dfrac{l+r}{2}$

首先我们对于区间 $[l,mid]$ 和 $[mid,r]$ 分别使用辛普森公式，得到答案为 $Lans,Rans$

那么如果 $|Lans+Rans-ans| \le 15 \times eps$ ，$eps$ 是需要的精度，我们就认为这个精度是可以接受的，不需要继续递归，直接返回 $ans$ 即可

其实我也不知道为什么是15倍的eps

否则，分左右两边递归，把答案相加返回即可

代码长这样：

double asr(double l,double r,double eps,double ans) {
    double mid = (l + r) / 2;
    double Lv = simpson(l,mid),Rv = simpson(mid,r);
    if(fabs(Lv + Rv - ans) <= 15 * eps) return Lv + Rv + (Lv + Rv - ans) / 15;
    return asr(l,mid,eps / 2,Lv) + asr(mid,r,eps / 2,Rv);
}

需要注意的是，区间缩小为之前的一半之后， $eps$ 也需要 $\times \dfrac{1}{2}$

然后就做完了

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

直接上板子即可，没有任何区别

#include<bits/stdc++.h>
using namespace std;

double a,b,c,d,L,R;

inline double f(double x) {return (c * x + d) / (a * x + b);}

inline double simpson(double l,double r) {
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}

double asr(double l,double r,double eps,double ans) {
    double mid = (l + r) / 2;
    double Lv = simpson(l,mid),Rv = simpson(mid,r);
    if(fabs(Lv + Rv - ans) <= 15 * eps) return Lv + Rv + (Lv + Rv - ans) / 15;
    return asr(l,mid,eps / 2,Lv) + asr(mid,r,eps / 2,Rv);
}

inline double ASR(double l,double r,double eps) {return asr(l,r,eps,simpson(l,r));}

int main () {
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
    printf("%.6lf\n",ASR(L,R,1e-6));
    return 0;
}

板子 P4526 【模板】自适应辛普森法2

这个稍微有一点点地方需要想想

它要我们求

$$\int_0^\infty x^{\frac{a}{x}-x}\,dx,a\le50$$

如果积分发散，输出 $orz$，保留 $5$ 位小数

显然当 $x = 0$ 的时候这东西是没意义的，所以变成积 $(0,+\infty)$

然后我们把这个函数丢到画图软件里面很快可以发现：

$a < 0$ 的时积分发散，$\lim_{x \to 0} f(x) \to +\infty$

同时我们可以发现，当 $a > 0$ 的时候，当 $x > 10$ 时函数就趋近于 $0$ 了，所以我们就可以把积分的区间变成 $(0,10]$

以上证明待补

它让我们保留 $5$ 位小数，当我的 $eps = 10^{-5}$ 时， $WA$ 了两个点，记得把 $eps$ 开得更小些

保险起见，我的右端设的是 $12$，$eps=10^{-8}$，注意左端不能取到 $0$ ！

#include<bits/stdc++.h>
using namespace std;

double a;

inline double f(double x) {return pow(x,(a / x) - x);}

inline double simpson(double l,double r) {
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}

double asr(double l,double r,double eps,double ans) {
    double mid = (l + r) / 2;
    double Lv = simpson(l,mid),Rv = simpson(mid,r);
    if(fabs(Lv + Rv - ans) <= 15 * eps) return Lv + Rv + (Lv + Rv - ans) / 15;
    return asr(l,mid,eps / 2,Lv) + asr(mid,r,eps / 2,Rv);
}

inline double ASR(double l,double r,double eps) {return asr(l,r,eps,simpson(l,r));}

int main () {
    scanf("%lf",&a);
    if(a < 0) return puts("orz"),0;
    printf("%.5lf\n",ASR(1e-8,12,1e-8));
    return 0;
}

例题

咕咕咕