Rings (AOJ2562)

クソコード博覧会.
3次元はだいたい球とかで線分とかとの話になって二次式を解くみたいなのになりがちだが、これはまあまあおもしろかった(設定が簡潔なのもよい)とりあえずきれいに回転させたあとz=0平面との交わりを出して円の内外にわかれてるか判定.交わりはとりあえず平面同士の交直線を出してうまい距離のところを見る.ベクトルっぽく考えると多分楽.回転する時正規直交基底というのを完全に忘れてて無駄な計算をしてしまった・・・(転置するだけでOK)

#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define rep1(i,n) for(int i=1;i<=(int)(n);i++)
#define all(c) c.begin(),c.end()
#define pb push_back
#define fs first
#define sc second
#define show(x) cout << #x << " = " << x << endl
#define chmin(x,y) x=min(x,y)
#define chmax(x,y) x=max(x,y)
using namespace std;
typedef double D;
typedef vector<D> vd;
typedef vector<vd> mat;
typedef complex<D> P;
typedef pair<P,P> L;
typedef vector<P> Pol;
struct C{P p;D r;};
D eps=1e-9;
bool eq(D a, D b) { return abs(a-b)<eps;}
inline int large(D a,D b){return eq(a,b)?0:(a>b?1:-1);}
D dot(P a,P b){return real(conj(a)*b);}
inline P perp(L l, P p){
	D t=dot(p-l.fs,l.fs-l.sc)/norm(l.fs-l.sc);
	return l.fs+t*(l.fs-l.sc);
}
inline D dLP(L l, P p) { return abs(perp(l,p)-p);}
inline Pol intCL(C c,L l){
	int x=large(dLP(l,c.p),c.r);
	Pol ret;
	if(x==1) return ret;
	if(x==0){
		ret.pb(perp(l,c.p));
		return ret;
	}
	D d=dLP(l,c.p);
	if(d<eps){
		ret.pb(c.p+(l.fs-l.sc)/abs(l.fs-l.sc)*c.r);
		ret.pb(c.p-(l.fs-l.sc)/abs(l.fs-l.sc)*c.r);
		return ret;
	}
	D theta=acos(d/c.r);
	P p=perp(l,c.p);
	ret.pb(c.p+(p-c.p)/d*polar(c.r,theta));
	ret.pb(c.p+(p-c.p)/d*polar(c.r,-theta));
	return ret;
}
vd pro(mat m,vd v){
	vd a(3,0);
	rep(i,3) rep(j,3) a[i]+=m[i][j]*v[j];
	return a;
}
int main(){
	D x,y,z,a,b,c,d,e,f,p,q,r,A,B,C_,D_,E,F;
	cin>>x>>y>>z>>a>>b>>c>>d>>e>>f;
	cin>>p>>q>>r>>A>>B>>C_>>D_>>E>>F;
	p-=x,q-=y,r-=z,x=0,y=0,z=0;
	D g=b*f-c*e,h=c*d-a*f,i=a*e-b*d;
	mat M={
		{e*i-f*h,g*f-d*i,d*h-g*e},
		{h*c-b*i,a*i-g*c,g*b-a*h},
		{b*f-e*c,d*c-a*f,a*e-b*d}
	};
	vd v=pro(M,vd{A,B,C_}),u=pro(M,vd{D_,E,F}),o=pro(M,vd{p,q,r});
	if(eq(u[2],0)) swap(v,u);
	if(eq(u[2],0)){
		puts("NO");
		return 0;
	}
	P p1=P(0,-o[2]/u[2]),p2=P(1,-(o[2]+v[2])/u[2]);
	L l=L(p1,p2);
	Pol pol=intCL(C{P(0,0),1},l);
	if(pol.size()<=1){
		puts("NO");
		return 0;
	}
	int cnt=0;
	for(P p:pol){
		D s=p.real(),t=p.imag();
		D xx=o[0]+s*v[0]+t*u[0],yy=o[1]+s*v[1]+t*u[1];
		if(xx*xx+yy*yy<1) cnt++;
	}
	if(cnt==1) puts("YES");
	else puts("NO");
}