被faebdc的shadow坑惨了。在考场上才发现自己的3D计算几何知识基本为0。

旋转

已知点 $P=(x,y,z)$ ，把它沿着转动 $\theta^{\circ}$ （从与向量相对的角度看过去，点是逆时针旋转的）首先可以给出绕三个轴的旋转矩阵

R_x(\theta)=\begin{bmatrix} 1&0&0\\ 0&\cos\theta&-\sin\theta\\ 0&\sin\theta&\cos\theta \end{bmatrix}

R_y(\theta)=\begin{bmatrix} \cos\theta&0&\sin\theta\\ 0&1&0\\ -\sin\theta&0&\cos\theta \end{bmatrix}

$R_y$ 有些不同，因为 $R_y$ 的变化中 $z$ 等价成了二维的 $x$ 轴， $x$ 被等价成了 $y$ 轴

R_z(\theta)=\begin{bmatrix} \cos\theta&-\sin\theta&0\\ \sin\theta&\cos\theta&0\\ 0&0&1 \end{bmatrix}

对于这三个特殊的旋转轴我们已经有了解决方案了，那么对于任意的轴要怎么办呢？我们可以将这个旋转分解：将整个坐标轴旋转，使得旋转轴和 $z$ 轴重合再将点 $P$ 绕 $z$ 轴旋转 $\theta$ 角再将整个坐标轴旋转回原位

算了直接上结论。。Miskcoo

\begin{bmatrix} \cos \theta + x^2(1 - \cos \theta) & xy(1 - \cos \theta) - z\sin \theta & xz(1 - \cos \theta) + y\sin \theta \\ yx(1 - \cos \theta) + z\sin \theta & \cos \theta + y^2(1 - \cos \theta) & yz(1 - \cos \theta) - x\sin \theta \\ zx(1 - \cos \theta) - y\sin \theta & zy(1 - \cos \theta) + x\sin \theta & \cos \theta + z^2(1 - \cos \theta) \end{bmatrix}

还是很优美的。

线面相交

设平面的法向量为，平面上有一点 $P$ 。有直线。

设交点为 $C$ ，则

解得

当的时候，显然交点不存在。

这种法线的操作对于二维的线线相交依然适用，比面积法或者体积法显然得多，好记得多，而且数值稳定性并没有损失。

一道题

附上一道练手题目：Shadow

这道题目存在更简洁的做法，但是为了强行套模板，我选择了一种非常难受的转化：把坐标系强行转化成地面是坐标平面，转化成把 $(0,0,0),B,C$ 构成的三角形转到 $xOy$ 面上。详见代码，实在是恶心得不想多说。

写完1A好评。

#include <cstdio>
struct Istream {
	template <class T> 
	Istream &operator >>(T &x) {
		static char ch;static bool neg;
		for(ch=neg=0;ch<'0' || '9'<ch;neg|=ch=='-',ch=getchar());
		for(x=0;'0'<=ch && ch<='9';(x*=10)+=ch-'0',ch=getchar());
		x=neg?-x:x;
		return *this;
	}
	Istream &operator >>(double &x) {
		scanf("%lf",&x);
		return *this;
	}
}is;
struct Ostream {
	template <class T>
	Ostream &operator <<(T x) {
		x<0 && (putchar('-'),x=-x);
		static char stack[233];static int top;
		for(top=0;x;stack[++top]=x%10+'0',x/=10);
		for(top==0 && (stack[top=1]='0');top;putchar(stack[top--]));
		return *this;
	}
	Ostream &operator <<(double const &x) {
		printf("%.2lf",x);
		return *this;
	}
	Ostream &operator <<(char ch) {
		putchar(ch);
		return *this;
	}
}os;
#include <cmath>
#include <cassert>
#include <algorithm>
const int MAXN=100+11;
const double eps=1e-7,pi=acos(-1);
inline int sgn(double const &x) { return (x>eps)-(x<-eps); }
#define SELF_OPER(Type,op) template <class T> Type &operator op##=(T const &x) { return *this=*this op x; } 
namespace _3D {
typedef struct Point Vector;
struct Point {
	double x,y,z;
	Point() {}
	Point(double const &x,double const &y,double const &z):x(x),y(y),z(z) {}
	double Length() const { return sqrt(x*x+y*y+z*z); }
	friend double Inner(Vector const &a,Vector const &b) { return a.x*b.x+a.y*b.y+a.z*b.z; }
	friend Vector Outer(Vector const &a,Vector const &b) { return Vector(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x); }
	friend Point operator +(Point const &lhs,Vector const &rhs) { return Point(lhs.x+rhs.x,lhs.y+rhs.y,lhs.z+rhs.z); }
	friend Vector operator -(Point const &lhs,Point const &rhs) { return Vector(lhs.x-rhs.x,lhs.y-rhs.y,lhs.z-rhs.z); }
	friend Vector operator *(Vector const &lhs,double const &rhs) {	return Vector(lhs.x*rhs,lhs.y*rhs,lhs.z*rhs); }
	friend Vector operator /(Vector const &lhs,double const &rhs) {	return Vector(lhs.x/rhs,lhs.y/rhs,lhs.z/rhs); }
	SELF_OPER(Point,+)
	SELF_OPER(Point,-)
	SELF_OPER(Point,*)
	SELF_OPER(Point,/)
	friend double Angle(Vector const &v1,Vector const &v2) { return asin(Outer(v1,v2).Length()/(v1.Length()*v2.Length())); }
	friend Istream &operator >>(Istream &is,Point &p) {
		is>>p.x>>p.y>>p.z;
		return is;
	}
	friend Ostream &operator <<(Ostream &os,Point const &p) {
		os<<'('<<p.x<<','<<p.y<<','<<p.z<<')';
		return os;
	}
}p[MAXN];
struct Line {
	Point p;
	Vector v;
	Line() {}
	Line(Point const &p,Vector const &v):p(p),v(v) {}
};
Point Rotate(Line const &ln,Point p,double const &rad) {
	static double co,si;
	static Vector a;
	si=sin(rad),co=cos(rad);
	a=ln.v/ln.v.Length();
	double T[][3]={{co+a.x*a.x*(1-co),a.x*a.y*(1-co)-a.z*si,a.x*a.z*(1-co)+a.y*si},
				  {a.y*a.x*(1-co)+a.z*si,co+a.y*a.y*(1-co),a.y*a.z*(1-co)-a.x*si},
				  {a.z*a.x*(1-co)-a.y*si,a.z*a.y*(1-co)+a.x*si,co+a.z*a.z*(1-co)}};
	p-=ln.p;
	p=Point(p.x*T[0][0]+p.y*T[1][0]+p.z*T[2][0],
			p.x*T[0][1]+p.y*T[1][1]+p.z*T[2][1],
			p.x*T[0][2]+p.y*T[1][2]+p.z*T[2][2]);
	return p+=ln.p;
}
struct Plane {
	Point p;Vector nv;
	Plane(Point const &p1,Point const &p2,Point const &p3):p(p1),nv(Outer(p2-p1,p3-p1)) {}
};
const Plane GROUND(Point(0,0,0),Point(1,0,0),Point(0,1,0));
Point Cross(Line ln,Plane pln) {
	double d1=-Inner(pln.nv,ln.p-pln.p),d2=Inner(pln.nv,ln.v);//d1/d2
	if(d2==0) {
		throw;
	} else {
		return ln.p+ln.v*(d1/d2);
	}
}
}
namespace _2D {
typedef struct Point Vector;
struct Point {
	double x,y;
	Point() {}
	Point(double const &x,double const &y):x(x),y(y) {}
	friend bool operator <(Point const &lhs,Point const &rhs) {	return lhs.x==rhs.x?lhs.y<rhs.y:lhs.x<rhs.x; }
	friend Vector operator +(Point const &lhs,Vector const &rhs) { return Vector(lhs.x+rhs.x,lhs.y+rhs.y); }
	friend Vector operator -(Point const &lhs,Point const &rhs) { return Vector(lhs.x-rhs.x,lhs.y-rhs.y); }
	friend Vector operator *(Vector const &lhs,double const &rhs) { return Vector(lhs.x*rhs,lhs.y*rhs); }
	friend Vector operator /(Vector const &lhs,double const &rhs) {	return Vector(lhs.x/rhs,lhs.y/rhs); }
	friend double Inner(Vector const &a,Vector const &b) { return a.x*b.x+a.y*b.y; }
	friend double Outer(Vector const &a,Vector const &b) { return a.x*b.y-a.y*b.x; }
}p[MAXN];
struct Line {
	Point p;Vector v;
	Line() {}
	Line(Point const &p,Vector const &v):p(p),v(v) {}
};
Point Cross(Line const &ln1,Line const &ln2) {
	Vector norm(-ln2.v.y,ln2.v.x);
	double d1=-Inner(norm,ln1.p-ln2.p),d2=Inner(ln1.v,norm);
	return ln1.p+ln1.v*(d1/d2);
}
double HullSize(Point p[],int n) {
	std::sort(p+1,p+1+n);
	static Point stack[MAXN];
	int top=0;
	for(int i=1;i<=n;++i) {
		while(top>=2 && Outer(stack[top]-stack[top-1],p[i]-stack[top])<=0) {
			top--;
		}
		stack[++top]=p[i];
	}
	int mid=top;
	for(int i=n;i>=1;--i) {
		while(top>=mid+1 && Outer(stack[top]-stack[top-1],p[i]-stack[top])<=0) {
			top--;
		}
		stack[++top]=p[i];
	}
	if(top) {
		top--;
	}
	double res=0;
	for(int i=1;i<=top;++i) {
		res+=Outer(stack[i],stack[i==top?1:i+1]);
	}
	return res/2;
}
}
using namespace _3D;
void Transform(Point g[],Point &light,Point p[],int n) {
	Point &o=g[0];
	for(int i=1;i<=n;++i) {
		p[i]-=o;
	}
	light-=o;
	g[2]-=o;
	g[1]-=o;
	g[0]-=o;
	Line axis;double angle;
	axis=Line(o,Vector(-g[1].y,g[1].x,0));
	angle=sgn(g[1].x)==0 && sgn(g[1].y)==0?pi/2:Angle(g[1],Vector(g[1].x,g[1].y,0));
	if(sgn(Rotate(axis,g[1],angle).z)!=0) {
		angle*=-1;
	}
	for(int i=1;i<=n;++i) {
		p[i]=Rotate(axis,p[i],angle);
	}
	light=Rotate(axis,light,angle);
	g[1]=Rotate(axis,g[1],angle);
	g[2]=Rotate(axis,g[2],angle);
	assert(sgn(g[1].z)==0);

	_2D::Point _crs=_2D::Cross(_2D::Line(_2D::Point(g[2].x,g[2].y),_2D::Vector(-g[1].y,g[1].x)),_2D::Line(_2D::Point(0,0),_2D::Point(g[1].x,g[1].y)));
	Point crs=Point(_crs.x,_crs.y,0);
	axis=Line(o,g[1]);
	angle=sgn(g[2].x-crs.x)==0 && sgn(g[2].y-crs.y)==0?pi/2:Angle(g[2]-crs,Point(g[2].x,g[2].y,0)-crs);
	if(sgn(Rotate(axis,g[2],angle).z)!=0) {
		angle*=-1;
	}
	for(int i=1;i<=n;++i) {
		p[i]=Rotate(axis,p[i],angle);
	}
	light=Rotate(axis,light,angle);
	g[2]=Rotate(axis,g[2],angle);
	assert(sgn(g[2].z)==0);
}

int main() {
#ifdef DEBUG
	freopen("input","r",stdin);
#else
	freopen("shadow.in","r",stdin);
	freopen("shadow.out","w",stdout);
#endif
	Point g[3];is>>g[0]>>g[1]>>g[2];
	Point light;is>>light;
	int n;is>>n;
	for(int i=1;i<=n;++i) {
		is>>p[i];
	}
	Transform(g,light,p,n);
	Plane pln(g[0],g[1],g[2]);
	for(int i=1;i<=n;++i) {
		Point crs=Cross(Line(light,p[i]-light),pln);
		_2D::p[i]=_2D::Point(crs.x,crs.y);
	}
	os<<_2D::HullSize(_2D::p,n)<<'\n';
	return 0;
}