NEERC 2008 - A - Aerodynamics

#include <bits/stdc++.h>

using namespace std;

#define mk make_pair
#define pb push_back
#define fi first
#define se second

typedef pair<int, int> ii;
typedef long long ll;
const double EPS = 1e-9;
const double PI = acos(-1.0);

ll gcd (ll a, ll b) {
	if (!b)
		return a;
	else
		return gcd(b, a%b);
}

class Point {
public:
	double x, y;

	Point () { }

	Point (double x, double y) {
		this->x = x;
		this->y = y;
	}

	/**/
	bool operator == (const Point &b) const {
		return (abs (x - b.x) < EPS and abs (y - b.y) < EPS);
	}

	/**/
	bool operator < (const Point &b) const {
		return ((x < b.x) or ((x == b.x) and y < b.y));
	}

	//Produto vetorial
	double operator ^ (const Point &b) const {
		return (this->x * b.y) - (this->y * b.x); 
	}

	//Produto escalar
	double operator * (const Point &b) const {
		return (this->x * b.x) + (this->y * b.y);
	}

	/**/
	Point operator * (double k) {
		return Point (x*k, y*k);
	}

	Point operator / (double k) {
		if (k == 0.0) cout << "Class Point (operador /): divisao por zero" << endl;
		return Point (x/k, y/k);
	}
	
	/**/
	Point operator + (const Point &b) const {
		return Point (x + b.x, y + b.y);
	}

	/**/
	Point operator - (const Point &b) const {
		return Point (x - b.x, y - b.y);
	}

	/**/
	double len () {
		return sqrt (x*x + y*y);
	}

	/*Distancia ponto a ponto*/
	double dpp (const Point &b) {
		return ((*this)-b).len();
	}

	/*Projecao*/
	double proj (Point &b) {
		return ((*this)*b)/(b.len()*b.len());
	}

	Point norm () {
		return Point (x/this->len(), y/this->len());
	}

	/*Retorna o vetor perpendicular ao vetor (0,0) -> (Point)
	Sentido clockwise*/
	Point perp () {
		return Point (this->y, -1.0 * this->x);
	}

	// Distancia do ponto p ao segmento ab, tambem retorna por 
	// referencia o ponto (c) do segmento mais perto de p
	double distToLine (const Point a, const Point b, Point& c) {
		// formula: c = a + u * ab
		Point p = *this;
		if (a == b) return p.dpp(a);
		Point ap = Point(p - a), ab = Point(b - a);
		double u = ap.proj(ab);
		if (u < 0.0) u = 0.0;
		if (u > 1.0) u = 1.0;
		c = a + ab * u;
		return p.dpp(c);
	}

	/**/
	Point rotaciona (double ang) {
		double c = cos(ang), s = sin(ang);
		double X = x*c - y*s;
		double Y = x*s + y*c;
		return Point(X,Y);
	}

	/*area de um poligono concavo ou convexo
	dado vetor de vertices ordenados clockwise ou
	counter clockwise*/
	static double area (vector <Point> v) {
		double area = 0.0;
		for (int i = 0; i < (int)v.size(); i++) 
			area += v[i] ^ v[(i+1)%v.size()];
		return abs(area/2.0);
	}
	
	/*Angulo em forma de fracao reduzida entre o vetor Op (p é o ponto)
	 e o eixo x, se paralelo ao eixo x retorna (1,0) ou (-1,0)
	 se paralelo ao eixo y retorna (0,1) ou (0,-1)
	 SÓ FUNCIONA PARA PONTOS INTEIROS*/ 
	static ii ang (Point p) {
		int a = p.x, b = p.y;
		if (a == 0) return mk(0, b/abs(b));
		else if (b == 0) return mk(a/abs(a), 0);
		return mk(a/gcd(a,b), b/gcd(a,b));
	}

	/* return counter clock points of convex hull*/
	static vector <Point> convex_hull (vector <Point> p) {
		if (p.size() <= 2) return p;

		int n = p.size(), k = 0;
		vector <Point> H(2*n);

		sort(p.begin(), p.end());

		for (int i = 0; i < n; i++) {
			while (k >= 2 and ((H[k-1]-H[k-2])^(p[i]-H[k-2])) <= 0.0) k--;
			H[k++] = p[i];
		}

		for (int i = n-2, t = k+1; i >= 0; i--) {
			while (k >= t and ((H[k-1]-H[k-2])^(p[i]-H[k-2])) <= 0.0) k--;
			H[k++] = p[i];
		}

		H.resize(k-1);

		return H;
	}
	
	/*Interseccao de dois vetores somandos a pontos*/
	static double inter (Point p1, Point v1, Point p2, Point v2) {
		if (abs(v2 ^ v1) >= EPS) {
			Point c = p1 - p2;
			return (c ^ v2)/(v2 ^v1);
		} else {
			cout << "Class Point (funcao inter): retas paralelas" << endl;
			cout << "Talvez deva ajustar o EPS" << endl;
		}
		return 0.0;
	}

	/*Retorna o retangulo (pontos em anti clockwise) que tem a menor valor
	min(Xmax -Xmin, Ymax - Ymin)*/
	static vector <Point> minRetangulo (vector <Point> v) {
		vector <Point> at;
		at.pb(Point(-1e18, -1e18));
		at.pb(Point(1e18, -1e18));
		at.pb(Point(1e18, 1e18));
		at.pb(Point(-1e18, 1e18));
		v = convex_hull(v);
		int n = v.size();
		for (int i = 0; i < n; i++) {
			int j = (i+1)%n;
			Point vec = v[j] - v[i];
			double ang = atan2(vec.y, vec.x);
			vector <Point> ve;
			for (int j = 0; j < n; j++) {
				ve.pb(v[j].rotaciona(ang));
			}
			double minx = DBL_MAX, miny = DBL_MAX, maxx = -DBL_MAX, maxy = -DBL_MAX;
			for (int j = 0; j < n; j++) {
				if (ve[j].x < minx) minx = ve[j].x;
				if (ve[j].x > maxx) maxx = ve[j].x;
				if (ve[j].y < miny) miny = ve[j].y;
				if (ve[j].y > maxy) maxy = ve[j].y;
			}
			double mini = min(maxx - minx, maxy - miny);
			if (mini < min(at[2].x - at[0].x, at[2].y - at[0].y)) {
				at.clear();
				at.pb(Point(minx, miny));
				at.pb(Point(maxx, miny));
				at.pb(Point(maxx, maxy));
				at.pb(Point(minx, maxy));
			}
		}

		return at;
	}
};

ostream &operator<<(ostream &os, Point const &p) {
	return os << p.x << " " << p.y;
}

class Circle {
public:
	Point c;
	double r;

	Circle () {}

	Circle (Point c, double r) : c(c), r(r) {}

	/*Interseccao de dois circulos
	OBS: se ha infinitas interseccoes retorna o vetor vazio
	OBS: se existe só um ponto retorna 2 pontos iguais*/
	vector <Point> intersect (Circle b) {
		vector <Point> ret;
		Point c1 = this->c, c2 = b.c;
		double r1 = this->r, r2 = b.r;

		if (c1 == c2) return ret;

		double d = c1.dpp(c2);
		
		if (d > r1 + r2 + EPS) return ret;
		if (d + EPS < abs(r1 - r2)) return ret;

		double a = (r1*r1 - r2*r2 + d*d)/(2.0*d);
		double h = sqrt(max(0.0, r1*r1 - a*a));

		Point pc = c1 + ((c2 - c1)*a)/d;

		/*X EH MENOS E Y EH MAIS*/
		double x = pc.x - ((h*(c2.y - c1.y))/d);
		double y = pc.y + ((h*(c2.x - c1.x))/d);
		ret.pb(Point(x,y));

		x = pc.x + ((h*(c2.y - c1.y))/d);
		y = pc.y - ((h*(c2.x - c1.x))/d);
		ret.pb(Point(x,y));

		return ret;
	}

};

/* Get area of a nondegenerate triangle, with sides a, b, c
*/
double getAreaTriangle (double a, double b, double c) {
	double p = (a + b + c)/2.0;
	return sqrt (p * (p - a) * (p - b) * (p - c));
}



/* return counter clock points of convex hull*/
static vector <Point> convex_hull (vector <Point> p) {
	if (p.size() <= 2) return p;

	int n = p.size(), k = 0;
	vector <Point> H(2*n);

	sort(p.begin(), p.end());

	for (int i = 0; i < n; i++) {
		while (k >= 2 and ((H[k-1]-H[k-2])^(p[i]-H[k-2])) <= 0.0) k--;
		H[k++] = p[i];
	}

	for (int i = n-2, t = k+1; i >= 0; i--) {
		while (k >= t and ((H[k-1]-H[k-2])^(p[i]-H[k-2])) <= 0.0) k--;
		H[k++] = p[i];
	}

	H.resize(k-1);

	return H;
}

typedef struct Point3d {
	double x,y,z;
}Point3d;

bool doubleEquals(double a, double b){
	if(max(a-b, b-a) < 1e-6)
		return true;
	return false;
}


vector<Point3d> pts2;
vector<Point> planePts;
vector<Point> hull;

int main (void) {
	freopen("aerodynamics.in", "r", stdin);
	freopen("aerodynamics.out", "w", stdout);
	
	ios_base::sync_with_stdio(false);
	cout << fixed << setprecision(10);
	
	int n,zmin,zmax;
	cin >> n >> zmin >> zmax;
	
	vector<Point3d> pts;
	while(n--){		
		Point3d p;		
		cin >> p.x >> p.y >> p.z;
		pts.pb(p);
	}
	
	//cout << "zcut" << endl;
	
	for(int zCut = zmin; zCut <= zmax; zCut++){
		pts2.clear();
		hull.clear();
		planePts.clear();
		
		//cout << "3d pts" << endl;
		for(Point3d p : pts){
			Point3d pn;
			pn.x = p.x;
			pn.y = p.y;
			pn.z = p.z - zCut;
			pts2.pb(pn);
		}
		
		
		//for(Point3d p : pts2)
		//	cout << p.x << " " << p.y << " " << p.z << endl;
		
		//cout << "Plane pts" << endl;
		
		for(int i = 0; i < pts2.size(); i++){
			for(int j = i+1; j < pts2.size(); j++){
				Point3d p1 = pts2[i];
				Point3d p2 = pts2[j];
				
				//cout << "Point pair1: " << p1.x << " " << p1.y << " " << p1.z << endl;
				//cout << "Point pair2: " << p2.x << " " << p2.y << " " << p2.z << endl;
				if(!doubleEquals(p2.z, p1.z)){
					double alpha = p1.z/(p1.z-p2.z);
					
					//cout << "ALPHA: " << alpha << endl;
					
					if(alpha >= 0 && alpha <= 1){
						Point pt;
						pt.x = p1.x + alpha*(p2.x-p1.x);
						pt.y = p1.y + alpha*(p2.y-p1.y);
						planePts.pb(pt);
					}
					
				}
				else if(doubleEquals(p2.z, 0)){
					Point pt;
					pt.x = p1.x;
					pt.y = p1.y;
					planePts.pb(pt);
					
					Point pt2;
					pt2.x = p2.x;
					pt2.y = p2.y;
					planePts.pb(pt2);
				}
				else {
					//cout << "INVALID" << endl;
				}
			}
		}
		
		//for(Point p : planePts)
		//	cout << p.x << " " << p.y << endl;
		
		//cout << "Calling hull" << endl;
		
		if(planePts.size() > 2){
			hull = convex_hull(planePts);
				
			//for(Point p : hull)
			//	cout << p.x << " " << p.y << endl;
			
			double area = 0;
			for(int i = 0; i < hull.size()-1; i++){
				area += hull[i].x*hull[i+1].y - hull[i].y*hull[i+1].x;
			}
			area += hull[hull.size()-1].x*hull[0].y - hull[hull.size()-1].y*hull[0].x;
			area /= 2;
			
			cout << area << endl;
		}
		else
			cout << 0 << endl;
	}
	
	return 0;
}