前言

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

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

正题

参考博客

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

$${\int }_l^{r}f(x)dx$$

我们显然有两种方法：

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

$${\int }_l^{r}f(x)dx=\underset{n\to \mathrm{\infty }}{lim}\sum _{i=1}^n{\displaystyle \frac{r-l}{n}}\times f(l+{\displaystyle \frac{r-l}{n}}\times i)$$
2、牛顿-莱布尼茨公式法，令 $F^{\prime }(x)=f(x)$，那么显然有：

$${\int }_l^{r}f(x)dx=F(r)-F(l)$$

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

辛普森公式

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

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

那我们有这样一个东西：

${\int }_l^{r}f(x)dx\approx {\int }_l^{r}g(x)dx$

$={\int }_l^r(a{x}^2+bx+c)dx$

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

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

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

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

$={\displaystyle \frac{(r-l)}{6}}\times (2a{r}^2+2alr+2a{l}^2+3br+3bl+6c)$

$={\displaystyle \frac{(r-l)}{6}}\times (a{l}^2+bl+c+a{r}^2+br+c+a{r}^2+2alr+a{l}^2+2br+2bl+4c)$

$={\displaystyle \frac{(r-l)}{6}}\times (g(l)+g(r)+a{r}^2+2alr+a{l}^2+2br+2bl+4c)$

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

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

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

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

我们就得出了这个公式：

$${\int }_l^{r}f(x)dx\approx {\displaystyle \frac{(r-l)(f(l)+f(r)+4f({\displaystyle \frac{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={\displaystyle \frac{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 {\displaystyle \frac{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^{\mathrm{\infty }}{x}^{\frac{a}{x}-x}dx,a\le 50$$

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

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

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

$a<0$ 的时积分发散，$\underset{x\to 0}{lim}f(x)\to +\mathrm{\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;
}

例题

咕咕咕