自适应辛普森求助
查看原帖
自适应辛普森求助
631787
Watanabe楼主2023/12/18 11:51

rt,调了两年半的参了,是精度问题还是啥求教

#include<iostream>
#include<cmath>
#include<cstring>
#include<vector>
#include<algorithm>
#include<cstdio>
#define re register
#define ld long double
#define PDD pair<long double,long double> 
#define x first
#define y second

using namespace std;
const int N=1e3+10;
const ld eps=1e-12,INF=1e6+10;
int n;
struct Trh{
	PDD a,b,c;
}P[N];

inline int sign(ld x)
{
	if(fabs(x)<eps) return 0;
	if(x<0) return -1;
	return 1;
}
inline ld get(PDD p,PDD q,ld x0)
{
	if(!sign(p.x-q.x)) return INF;
	ld lk=(p.y-q.y)/(p.x-q.x);
	ld lb=p.y-p.x*lk;
	ld A=x0*lk+lb;
	if(p.x>q.x) swap(p,q);
	if(sign(x0-p.x)<0||sign(x0-q.x)>0) return INF;
	return A;
}
PDD q[N];
inline ld f(ld x)
{
	int cnt=0;
	for(re int i=1;i<=n;++i)
	{
		auto A=get(P[i].a,P[i].b,x);
		auto B=get(P[i].a,P[i].c,x);
		auto C=get(P[i].b,P[i].c,x);
		if(!sign(A-INF)&&!sign(B-INF)&&!sign(C-INF)) continue ;
		int tot=0;
		if(sign(A-INF)!=0) q[++cnt].x=A,++tot;
		if(sign(B-INF)!=0)
		{
			if(!tot) q[++cnt].x=B;
			else q[cnt].y=B;
		}
		if(sign(C-INF)!=0) q[cnt].y=C;
	}
	if(!cnt) return 0;
	ld st,ed;
	st=q[1].x,ed=q[1].y;
	ld res=0;
	for(re int i=2;i<=cnt;++i)
	{
		if(q[i].x<=ed) ed=max(ed,q[i].y);
		else
		{
			res+=ed-st;
			st=q[i].x,ed=q[i].y;
		}
	}
	res+=ed-st;
	return res;
} 
inline ld S(ld l,ld r)
{
	return fabs((r-l)*(f(l)+f((l+r)/2.0)*4.0+f(r))/6.0);
}
ld simpson(ld l,ld r,ld s,int KK)
{
	auto mid=(l+r)/2.0;
	auto sl=S(l,mid),sr=S(mid+eps,r);
	if(!sign(s-sl-sr)||KK>=5) return sl+sr;
	return simpson(l,mid,sl,KK+1)+simpson(mid+eps,r,sr,KK+1);
}
vector <ld> ax;
signed main()
{
	cin>>n;
	for(re int i=1;i<=n;++i) cin>>P[i].a.x>>P[i].a.y>>P[i].b.x>>P[i].b.y>>P[i].c.x>>P[i].c.y;
	for(re int i=1;i<=n;++i) ax.push_back(P[i].a.x),ax.push_back(P[i].b.x),ax.push_back(P[i].c.x);
	sort(ax.begin(),ax.end());
	ld ans=0;
	for(re int i=0;i+1<ax.size();++i)
	{
		if(!sign(ax[i]-ax[i+1])) continue ;
		ans+=simpson(ax[i],ax[i+1],S(ax[i],ax[i+1]),1);
	}
	printf("%.2Lf",ans);
	return 0;
}
2023/12/18 11:51
加载中...