自适应辛普森了解一下

A 了这道超短的紫题，来发表一下本身的一些想法...

简单介绍辛普森这玩意儿

不如先学学泰勒展开？

首先泰勒展开你们都据说过吧？【雾

没据说过？安利某知乎回答：苍老师教你如何更好地记忆泰勒展开

而后你就知道了，泰勒展开实际上是对于某个函数在一个点不断去高阶求导，而后用求导获得的信息构造一个多项式，使得这个多项式在必定范围内几乎和原函数拟合（能够理解为接近重合的意思吧...）

类比到辛普森？

那么其实自适应辛普森也是相似的道理，只不过它是用了分治的方法去构造这个 假 拟合 多项式（其实就是二次函数，究其缘由应该就是这玩意儿比较好积分而且存在弧度，容易拟合吧...）

自适应是个什么鬼？

至于为何前面加了个自适应呢？由于咱们考虑分治是要有终止条件的（像对于一个序列分治的话就是分治的区间长度为 1 时中止），可是咱们这里是在实数域上拟合一个多项式啊，不存在什么规定的终止点...

因而咱们考虑怎样去设置终止条件？那固然是给定一个偏差范围（好比 1e-6），而后若是拟合函数和原函数大多数点的 y 值偏差不超过这个范围就终止分治，直接拿当前的函数去积分就行了，固然具体怎么比较的先别管（你看到下面以后会发现根本不须要比较两个函数 2333 ）

而后咱们发现这个偏差越大答案越不许确，越小跑的越慢，那么咱们就要调整这个偏差范围，而后这样的过程就有点像"自适应"了

拟合函数怎么构造？

咱们先将原函数约等于成一个二次多项式：

\[f(x)≈Ax^2+Bx+C\]

而后题目要咱们求的东西也转化一下：

\[ANS=\int_a^b f(x)dx\]

\[≈\int_a^b Ax^2+Bx+C\]

\[={A\over 3}(b^3-a^3) + {B\over 2} (b^2-a^2) +C(b-a)\]

\[={(b-a)\over 6} [ 2A(b^2+ab+a^2) + 3B(b+a)+6C]\]

\[={(b-a)\over 6} [ (Aa^2 + Ba+C) + (Ab^2+Bb+C)+4(~A~({a+b\over 2})^2 + B ({a+b\over2})+C)]\]

\[≈{(b-a)\over 6} [~ f(a) + f(b) + 4~f({a+b\over2})~]\]

恩？你问我推导哪里来的？ 大佬%%%

咱们发现最后约回去了，A B C 都不见了，并且这时候式子里面已经没有积分了...

因此怎么构造 A、B、C 仍是没讲？

没错，咱们发现上面其实根本不须要用到这个拟合二次函数的具体系数 A、B、C ，只须要将 f 的式子带入计算求值就行了

可是这样的话咱们发现以前说的终止条件不见了，由于咱们的拟合函数已经不须要了（并且找出来也很麻烦...）

其实咱们只须要比较当前区间 \([a,b]\) 不分治时候的答案 \(ans\) 和 分治下的答案 \(ansL+ansR\) 的偏差是否超过了咱们设置的精度范围就行了

至于正确性？（我怎么知道）

由于咱们知道这个分治确定是层数越多越精确的，因此分治的结果精确度确定高于当前 ans 的精确度，那么咱们卡卡精度就能让答案在偏差容许范围内了

FAQ: mmp 原来真的不用比较原函数和拟合函数...

若是要比较的话又能怎么比较呢？反正我是不会的...

code

代码超级短！相应的咱们能够知道这里咱们只须要将 f 函数改一改就能解决其余函数的积分问题了（必定精度下）

//by Judge
#include<bits/stdc++.h>
#define db double
using namespace std;
db a,b,c,d,l,r;
inline db f(db x){return (c*x+d)/(a*x+b);}
inline db simpson(db l,db r){ db mid=(l+r)/2;
    return (f(l)+f(r)+4*f(mid))*(r-l)/6;
}
db asr(db l,db r,db eps,db ans){
    db mid=(l+r)/2,L=simpson(l,mid),R=simpson(mid,r);
    if(fabs(L+R-ans)<=eps) return L+R+(L+R-ans);
    return asr(l,mid,eps/2,L)+asr(mid,r,eps/2,R);
}
inline db asr(db l,db r,db eps){
    return asr(l,r,eps,simpson(l,r));
}
int main(){ scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    return !printf("%.6lf\n",asr(l,r,1e-8));
}