自适应Simpson积分求助
查看原帖
自适应Simpson积分求助
285617
黑影洞人楼主2022/3/29 17:06
#include<cstdio>
#include<algorithm>
#include<vector>
#include<stack>
#include<cmath>
#define N 114
#define d double
#define pii pair<double ,double >
#define h first
#define s second
#define suky {114514,1919810}
using namespace std;
const d eps=1e-9;
const d inf=1e7;
int n,m;
d lll=inf,rrr=-inf;
bool mxe(d a,d b){
	return (a-b)>eps;
}
d gy(d x,d y,d x2,d y2,d t){
	d l=min(x,x2),r=max(x,x2),px=0,pn=0;
	if(mxe(l,t)||mxe(t,r))return 114514;
	d k=(x-x2)/(y-y2),b=y-((x2-x)*x/(y2-y));//一次函数 
	return k*t+b;
}
struct tri{
	d x[4]={0},y[4]={0};
	d l=inf,r=-inf;
	int lp=0,rp=0;
	void csh(){
		for(int i=1;i<=3;i++){
			if(l>x[i])lp=i,l=x[i];
			if(r<x[i])rp=i,r=x[i];
		}
	 }
	 pii up(d t){
		if(mxe(l,t)||mxe(t,r))return suky;
		d ya=gy(x[1],y[1],x[2],y[2],t);
		d yb=gy(x[3],y[3],x[2],y[2],t);
		d yc=gy(x[1],y[1],x[3],y[3],t);
		if((ya==114514&&yb==114514)||(yb==114514&&yc==114514)||(ya==114514&&yc==114514))return suky;
		if(ya==114514)return {max(yb,yc),min(yb,yc)};
		if(yb==114514)return {max(ya,yc),min(ya,yc)};
		if(yc==114514)return {max(yb,ya),min(yb,ya)};
	}
}tri[N]; 
d f(d x){
	d ans=0;
	vector<pii> v;
	for(int i=1;i<=n;i++){
		pii tmp=tri[i].up(x);
		if(tmp.h==114514&&tmp.s==1919810)continue;
		v.push_back(tmp);
	}
	sort(v.begin(),v.end());
	stack<d> s;
	int p=0;
	for(int i=0;i<v.size();i++){
		s.push(v[i].s);
		while(v[p].s<v[i].h)p++;
		ans+=v[p].h-v[i].s;
		i=p;
	}
	return ans;
}
d simpson(d l,d r){
	return (f(l)+4*f((l+r)/2)+f(r))*(r-l)/6;
}
double calc(double l,double r,double sim){
	double mid=(l+r)/2;
	double ll=simpson(l,mid),rr=simpson(mid,r);
	if(fabs(ll+rr-sim)<=15*eps)return ll+rr+(ll+rr-sim)/15.0;
	return calc(l,mid,ll)+calc(mid,r,rr);
}
signed main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++){
		scanf("%lf%lf%lf%lf%lf%lf",&tri[i].x[1],&tri[i].y[1],&tri[i].x[2],&tri[i].y[2],&tri[i].x[3],&tri[i].y[3]);
		tri[i].csh();
		lll=min(lll,tri[i].l);
		rrr=max(rrr,tri[i].r);
	}
	printf("%lf\n",calc(lll,rrr,simpson(lll,rrr)));
	return 0;
}



2022/3/29 17:06
加载中...