[USACO2010 Mar 【starc】]半平面交、~我写的是O(n^2)的

【题目链接】http://ace.delos.com/ioiupload?init=1&a=uNRSh9JKRlk
【题目大意】
给定n<=300个形如ax+by+cz>=0的不等式。(将<=0的乘一个 -1来改变不等号方向)
求在其满足的前提下，给出m个询问，
看能不能确定另外给出的式子ax+by+cz是>0还是<0还是无法确定。
它很猥琐的还有常量条件。
x>0
y>0
z>0
x<=100y
x<=100z
y<=100x
y<=100z
z<=100x
z<=100y
就这么多。
【算法分析】
由于x,y,z>0，所以我们将不等式两边同除以z，不用变号。
那么不等式变成：a(x/z)+b(x/z)+c>=0
我们令X=x/z,Y=x/z。
那么不等式变为aX+bY+c>=0。
然后？
半平面交！
就是先对给出的条件全部变为这种限制，然后求出其半平面交。
然后这个半平面交的算法是我自己YY出来的,O(n^2)。。。小有成就感。
当然，众神牛估计要笑而不语了。
好了，半平面交完以后呢，就只要对每个给出的限制转化为直线。然后看这个直线是“穿过”这个半平面交。
如果穿过，表示两边的取值都可能，所以不能确定他的不等号方向。
否则，看就看它在哪一边就行了。
【其它】
好有好多好多小细节。。。然后我写成的程序，精度要卡得刚好才能过= =
如果要学习半平面交的话。。。还是别看我这。因为写得实在太水，太长了。
我看过师傅的模板，不禁为其短小精悍而叹服。
但是毕竟是自己推得算法。。。再水也是自己的= =
【CODE】
/*
ID:jie22221
TASK:starc
LANG:C++
*/
#include #include #include #include #define N 355
#define INF 100
#define eps 1e-6
#define seps 1e-2
using namespace std;
struct Line{double a,b,c;}LL[N];
struct Point{double x,y;}p[N*N],tp[N*N];
int n,Ln,m,tn;

void init(){
int i,A1,B1,C1,A2,B2,C2;
char ch;
cin >> Ln >> m;
for (i=1;i<=Ln;i++){
cin >> ch >> A1 >> B1 >> C1 >> A2 >> B2 >> C2;
if (ch==’J’){
LL[i].a=A1-A2;
LL[i].b=B1-B2;
LL[i].c=C1-C2;
}else{
LL[i].a=A2-A1;
LL[i].b=B2-B1;
LL[i].c=C2-C1;
}
}
Ln++;
LL[Ln].a=-1;
LL[Ln].b=100;
LL[Ln].c=0;
Ln++;
LL[Ln].a=100;
LL[Ln].b=-1;
LL[Ln].c=0;
n=4;
p[1].x=seps; p[1].y=seps;
p[2].x=INF; p[2].y=seps;
p[3].x=INF; p[3].y=INF;
p[4].x=seps; p[4].y=INF;
}

inline double f(Line L,Point pp){
return L.a*pp.x+L.b*pp.y+L.c;
}
inline double cj(Point A,Point B){return A.x*B.y-A.y*B.x;}
inline void dec(int &x){x–; if (x<1) x=n;}
inline void inc(int &x){x++; if (x>n) x=1;}

Line Get_Line(Point A,Point B){
Line res;
res.a=A.y-B.y;
res.b=B.x-A.x;
res.c=cj(A,B);
return res;
}

Point Cross_Point(Line A,Line B){
Point res;
res.x=(B.b*A.c-A.b*B.c)/(B.a*A.b-A.a*B.b);
res.y=(B.a*A.c-A.a*B.c)/(A.a*B.b-B.a*A.b);
return res;
}

double fabs(double x){
if (x<0) return -x;
return x;
}

void Cut(Line L){
int i,j,k,st=1;
bool flag=false;
for (i=1;i<=n;i++)
if (f(L,p[i])<-eps) flag=true;
if (!flag){
for (i=1;i<=n;i++) tp[i]=p[i];
tn=n;
return;
}
while (f(L,p[st])<-eps) st++;
for (i=st;;inc(i))
if (f(L,p[i])<-eps) break;
for (j=st;;dec(j))
if (f(L,p[j])<-eps) break;
tn=0;
for (k=j%n+1;k!=i;inc(k))
tp[++tn]=p[k];
dec(i);
if (fabs(f(L,tp[tn]))>eps){
Line tmp=Get_Line(p[i],p[i%n+1]);
tp[++tn]=Cross_Point(tmp,L);
}
if (fabs(f(L,tp[ 1]))>eps){
Line tmp=Get_Line(p[j],p[j%n+1]);
tp[++tn]=Cross_Point(tmp,L);
}
}

void half_plane(){
int i,j;
for (i=1;i<=Ln;i++){
Cut(LL[i]);
for (j=1;j<=tn;j++) p[j]=tp[j];
n=tn;
}
}

void solve(){
int A1,B1,C1,A2,B2,C2,t1,t2,i,j;
for (i=1;i<=m;i++){
t1=t2=0;
cin >> A1 >> B1 >> C1 >> A2 >> B2 >> C2;
Line L;
L.a=A1-A2;
L.b=B1-B2;
L.c=C1-C2;
for (j=1;j<=n;j++){
if (f(L,p[j])<-1e-8) t1++;
if (f(L,p[j])> 1e-8) t2++;
}
if (t1==0 && t2==0 || t1+t2!=n) puts("U"); else
if (!t1) puts("J"); else
if (!t2) puts("B"); else
puts("U");
}
}

int main(){
freopen("starc.in","r",stdin);
freopen("starc.out","w",stdout);
ios::sync_with_stdio(false);
init();
half_plane();
solve();
return 0;
}