自适应Simpson积分

用$g(x)=A*x^2+B*x+C$来代替原函数$f(x)$，设$g(x)$原函数为$G(x)$，显然$G(x)=\frac{1}{3}*A*x^3+\frac{1}{2}*B*x^2+C*x$

$\int_a^b f(x) dx \approx \int_a^b g(x) dx=G(b)-G(a)$

代入化简得

$$\int_a^b f(x) dx \approx \frac{a-b}{6}*[g(a)+g(b)+4*g(\frac{a+b}{2})]$$

二分区间，然后计算区间积分，通过期望容差来控制二分

如果$|S(a,mid)+S(mid,b)-S(a,b)|<15*\varepsilon $，则中止二分，并且认为

$$S(a,b)=S(a,m)+S(m,b)+\frac{S(a,mid)+S(mid,b)-S(a,b)}{15}$$

否则继续二分

其中$mid=\frac{a+b}{2},\varepsilon为期望容差$

参考博客：https://www.cnblogs.com/chy-2003/p/11644112.html

例题：

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

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;

const double INF = 20;

double a, b, c, d, l, r;

inline double sqr(double x)
{
    return x * x;
}

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

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

double simpson(double l, double r, double eps)
{
    double mid = (l + r) / 2;
    double st = calc(l, r), sl = calc(l, mid), sr = calc(mid, r);
    if (fabs(sl + sr - st) <= 15 * eps) return sl + sr + (sl + sr - st) / 15;
    return simpson(l, mid, eps / 2) + simpson(mid, r, eps / 2);
}

int main()
{
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &l, &r);
    printf("%.6lf\n", simpson(l, r, 4e-7));
    return 0;
}

luogu 4526 【模板】自适应辛普森法2

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;

const double INF = 20;

double a;

inline double sqr(double x)
{
    return x * x;
}

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

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

double simpson(double l, double r, double eps)
{
    double mid = (l + r) / 2;
    double st = calc(l, r), sl = calc(l, mid), sr = calc(mid, r);
    if (fabs(sl + sr - st) <= 15 * eps) return sl + sr + (sl + sr - st) / 15;
    return simpson(l, mid, eps / 2) + simpson(mid, r, eps / 2);
}

int main()
{
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    scanf("%lf", &a);
    if (a < 0) printf("orz\n");
    else printf("%.5lf\n", simpson(4e-7, INF, 4e-7));
    return 0;
}

hdu 1724 Ellipse

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;

int T;
double a, b, l, r;

inline double sqr(double x)
{
    return x * x;
}

double f(double x)
{
    double t = sqr(b) - sqr(b) / sqr(a) * sqr(x);
    return 2 * sqrt(t);
}

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

double simpson(double l, double r, double eps)
{
    double mid = (l + r) / 2;
    double st = calc(l, r), sl = calc(l, mid), sr = calc(mid, r);
    if (fabs(sl + sr - st) <= 15 * eps) return sl + sr + (sl + sr - st) / 15;
    return simpson(l, mid, eps / 2) + simpson(mid, r, eps / 2);
}

int main()
{
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    scanf("%d", &T);
    while (T--) {
        scanf("%lf%lf%lf%lf", &a, &b, &l, &r);
        printf("%.3lf\n", simpson(l, r, 1e-6));
    }
    return 0;
}