感觉就是骗分工具，就像模拟退火一样

神仙的博客

P4526自适应辛普森法2

#include<bits/stdc++.h>
typedef double db;
using namespace std;
const db eps=1e-7;
db a;
db f(db x)
{
	return pow(x,a/x-x);
}
db simpson(db l,db r)
{
	return (f(l)+f(r)+4*(f((l+r)/2)))*(r-l)/6;
}
db integral(db l,db r,db eps,db a)
{
	db mid=(l+r)/2;
	db L=simpson(l,mid),R=simpson(mid,r);
	if(abs(L+R-a)<=eps*15) return L+R+(L+R-a)/15;
	return integral(l,mid,eps/2,L)+integral(mid,r,eps/2,R);
}
int main()
{
	scanf("%lf",&a);
	if(a<0) printf("orz\n");
	else printf("%.5lf\n",integral(eps,20,eps,simpson(eps,20)));
	return 0;
} 

也可以用来求圆的面积并(尽管不是正解)

#include<bits/stdc++.h>
typedef double db;
using namespace std;
const int N=1100;
const db blo=10;
struct seg
{
	db l,r;
	seg(db x=0,db y=0):l(x),r(y){}
};
int fa[N];
db mx[N],mn[N];
int n;
db x[N],y[N],r[N];
map<double,double> hav;
int find(int x)
{
	if(fa[x]==x) return x;
	return fa[x]=find(fa[x]);
}
void merge(int a,int b)
{
	a=find(a),b=find(b);
	fa[a]=b;
	mx[b]=max(mx[b],mx[a]);
	mn[b]=min(mn[b],mn[a]);
	return;
}
int cmp(seg a,seg b)
{
	return a.l<b.l;
}
db f(db X)
{
	vector<seg> a;
	for(int i=1;i<=n;i++)
	{
		db d=abs(x[i]-X);
		if(r[i]<d) continue;
		d=sqrt(r[i]*r[i]-d*d);
		a.push_back(seg(y[i]-d,y[i]+d));
	}
	sort(a.begin(),a.end(),cmp);
	int siz=a.size();
	if(siz==0) return 0;
	db sum=0,l=a[0].l,r=a[0].r;
	for(int i=1;i<siz;i++)
	{
		if(a[i].l<=r) r=max(a[i].r,r);
		else
		{
			sum+=r-l;
			l=a[i].l;
			r=a[i].r;
		}
	}
	sum+=r-l;
	return sum;
}
db simpson(db l,db r)
{
	return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}
db integral(db l,db r,db eps,db a)
{
	db mid=(l+r)/2;
	db L=simpson(l,mid),R=simpson(mid,r);
	if(abs(L+R-a)<=eps*15) return L+R+(L+R-a)/15;
	return integral(l,mid,eps/1.5,L)+integral(mid,r,eps/1.5,R);
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++) scanf("%lf%lf%lf",&x[i],&y[i],&r[i]);
	for(int i=1;i<=n;i++) fa[i]=i,mn[i]=x[i]-r[i],mx[i]=x[i]+r[i];
	for(int i=1;i<=n;i++)
	{
		for(int j=i+1;j<=n;j++)
		{
			if(abs(x[i]-x[j])<=r[i]+r[j])
			{
				merge(i,j);
			}
		}
	}
	db ans=0;
	double cnt=0;
	for(int i=1;i<=n;i++) if(fa[i]==i) cnt+=1;
	for(int i=1;i<=n;i++) if(fa[i]==i) ans+=integral(mn[i],mx[i],1e-3/cnt,simpson(mn[i],mx[i]));
	printf("%.3lf\n",ans);
	return 0;
}