自适应辛普森积分求助!!
查看原帖
自适应辛普森积分求助!!
519384
Link_Cut_Y楼主2022/4/17 12:13
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 110;
const double eps = 1e-6;
int n;
PDD q[N];

struct Triangle
{
	PDD p[3];
}tr[N];

double f(double x)
{
	int cnt = 0;
	bool have_intersection = false;
	
	for (int i = 0; i < n; i ++ )
	{
		for (int az = 0; az < 3; az ++ )
		{
		    int j = az;
			int k = (j + 1) % 3; // 当前线段是 tr[j] - tr[k]
			if (tr[i].p[k].x < tr[i].p[j].x) swap(j, k);
			if (x >= tr[i].p[j].x && x <= tr[i].p[k].x) // 判断不经过线段的情况
			{
    			have_intersection = true; 
    			double x1 = tr[i].p[j].x, y1 = tr[i].p[j].y;
    			double x2 = tr[i].p[k].x, y2 = tr[i].p[k].y;
    			
    			double k1 = (y2 - y1) / (x2 - x1);
    			double b = y1 - k1 * x1;
    			// 得到表达式 y = k1 * x + b
    			// 得到交点坐标 x0, y0 
    			double x0 = x, y0 = k1 * x0 + b;
    			if (q[cnt].x == 0) q[cnt].x = y0;
    			else q[cnt].y = y0;
			}
		}
		if (have_intersection) cnt ++ , have_intersection = false;
	}
	
	if (!cnt) return 0;
	for (int i = 0; i < cnt; i ++ )
	    if (q[i].x > q[i].y)
	        swap(q[i].x, q[i].y);
	        
	sort(q, q + cnt);
	
	double res = 0, st = q[0].x, ed = q[0].y;
	for (int i = 1; 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;
}

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

double query(double l, double r, double s)
{
	auto mid = (l + r) / 2;
	auto left = simpson(l, mid), right = simpson(mid, r);
	if (fabs(s - left - right) < eps) return left + right;
	return query(l, mid, left) + query(mid, r, right);
}

int main()
{
	scanf("%d", &n);
	
	double l = 1e8, r = -1e8;
	for (int i = 0; i < n; i ++ )
	{
		scanf("%lf%lf", &tr[i].p[0].x, &tr[i].p[0].y);
		scanf("%lf%lf", &tr[i].p[1].x, &tr[i].p[1].y);
		scanf("%lf%lf", &tr[i].p[2].x, &tr[i].p[2].y);
		
		for (int j = 0; j < 3; j ++ )
			l = min(l, tr[i].p[j].x),
			r = max(r, tr[i].p[j].x);
	}
	
	printf("%.2lf", query(l, r, simpson(l, r)));
	
	return 0;
}
2022/4/17 12:13
加载中...