Fork me on GitHub

Uoj#44.【清华集训2014】Breaking Bomber

Uoj#44. 【清华集训2014】Breaking Bomber

题意:

  • 给定元素矩阵$M$, $M$ 为一个$4\times N$的矩阵。
  • 给定一个四维列向量$x$ , 有两种要求:
    • 是否存在一个$N$维的列向量$w$, 满足存在$x=M\times w$。
    • 是否存在一个$N$维的列向量$w$, 满足存在 $x^T\times M\times w\ge 0$。
  • $N\le 5\times 10^4,M\le 10^5$。

题解:

  • 如果把每列抽象成一个四维空间的点,那么问题就是:
    • 给出一些四维空间中的点。
    • 询问一个点是否在这些点的凸包内。
    • 询问一个三维空间是否切到了这个凸包。
  • 看上去四维凸包不可做,但实际上可以转化为三维问题。
  • 对于一个点$(a,b,c,d)$变为$\frac{(a,b,c)}{a+b+c+d}$。
  • 对于一个三维空间$(A,B,C,D)$变为一个面$(A-D,B-D,C-D,D)$。
  • $Aa+Bb+Cc+Dd\ge 0=>(A-D)a+(B-D)b+(C-D)c+D(a+b+c+d)\\=>(A-D)\frac{a}{a+b+c+d}+(B-D)\frac{b}{a+b+c+d}+(C-D)\frac{c}{a+b+c+d}+D\ge 0$
  • 然后问题变为了:
    • 给出一些三维空间中的点。
    • 询问一个点是否在这些点的凸包内。
    • 询问一个面是否切到了这个凸包。
  • 听说有个结论:
    • 询问点是否在凸包中与询问半(超)平面是否切到半(超)平面交对偶。
    • 一个面$(A,B,C,D)$对偶成点$(\frac{A}{D},\frac{B}{D},\frac{C}{D})$。
  • 所以我们可以建出一开始所有点的凸包,然后将凸包上的面对偶成点,再求一遍凸包,询问一二就做法一样了。
  • 也就是说我们只需要求出三维凸包即可。
  • 然后有一种随机增量法,听说复杂度是$O(n\log n)$。
    • 随机加入点$p$。
    • 删掉所有可见面。
    • 加入$p$所在的新面。
    • 更新$p$可见面的信息和新面的可见点的信息。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
#include <bits/stdc++.h>
#define gc getchar()
#define sqr(x) ((x)*(x))
using namespace std;
const int N=50009;
const double eps=1e-9;
struct point
{
double x,y,z;
point(double x=0,double y=0,double z=0):x(x),y(y),z(z){}
double len()
{
return sqrt(x*x+y*y+z*z);
}
};
vector<point> points;
point operator +(const point &A,const point &B)
{
return point(A.x+B.x,A.y+B.y,A.z+B.z);
}
point operator -(const point &A,const point &B)
{
return point(A.x-B.x,A.y-B.y,A.z-B.z);
}
point operator *(const point &A,const double &B)
{
return point(A.x*B,A.y*B,A.z*B);
}
point operator *(const double &B,const point &A)
{
return point(A.x*B,A.y*B,A.z*B);
}
point operator /(const point &A,const double &B)
{
return point(A.x/B,A.y/B,A.z/B);
}
double operator *(const point &A,const point &B)
{
return A.x*B.x+A.y*B.y+A.z*B.z;
}
point cross(const point &A,const point &B)
{
return point(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);
}
point norm(point A)
{
return A/A.len();
}
struct face
{
point normal;
double d;
face(){}
face(point A,point B,point C)//normal*(x,y,z)=d
{
normal=norm(cross(B-A,C-A));
d=normal*A;//原点到平面距离
}
};
double face_point_dis(const face &A,const point &B)//normal.len=1
{
return A.normal*B-A.d;
}
bool see_face(const face &A,const point &B)//normal*point>d
{
return face_point_dis(A,B)>-eps;
}
point face_to_point(face A,point B)//将平面A关于B得到对偶点
{
A.d-=A.normal*B;
return A.normal/A.d;
}
struct Cpoint;
struct Cline;
struct Cface;
struct Cpoint
{
point p;
vector<Cface*> see,made;//see:看见的面,made:由其组成的面
};
struct Cline
{
Cface *f;//其所在的平面
Cpoint *p;//面上三点组成三角形中对应的点
Cline *back;//与其所在的平面以这条边相接的平面上的这条边
Cline *next_line();
};
struct Cface
{
face f;
Cpoint *killed;//被哪个点看见所删
vector<Cpoint*> seen;//可以看见这个平面的点集
Cline line[3];//平面的三条边
Cface(Cpoint *a,Cpoint *b,Cpoint *c):f(a->p,b->p,c->p),killed(NULL)//设置一个面的三条边的信息
{
line[0]=(Cline){this,b,NULL};
line[1]=(Cline){this,c,NULL};
line[2]=(Cline){this,a,NULL};
}
void set_line_back(Cline *a,Cline *b,Cline *c)//设置三条边的对应相接平面上的直线
{
line[0].back=a,a->back=line+0;
line[1].back=b,b->back=line+1;
line[2].back=c,c->back=line+2;
}
void find_see(const vector<Cpoint*> &now)
{
for (auto p:now)
if (p!=line[2].p&&see_face(f,p->p))//不是第一个点
seen.push_back(p),p->see.push_back(this);
}
void find_see(const vector<Cpoint*> &A,const vector<Cpoint*> &B)//从两个可能点集中选出可见点集
{
seen.resize(A.size()+B.size());
seen.clear();
__typeof(A.begin()) it1=A.begin();
__typeof(B.begin()) it2=B.begin();
Cpoint* Max=(Cpoint*)-1;
while (it1!=A.end()||it2!=B.end())
{
Cpoint* a=(it1==A.end())?Max:*it1;
Cpoint* b=(it2==B.end())?Max:*it2;
if (a<=b) it1++;
if (b<=a) it2++;
Cpoint* Min=min(a,b);
if (Min==line[2].p) continue;//不是第一个点
if (a==b||see_face(f,Min->p))
seen.push_back(Min),Min->see.push_back(this);
}
}
};
Cline* Cline::next_line()
{
return this->f->line+2==this?this->f->line:this+1;//这个面对应的下一条直线
}
struct Hull3D
{
vector<Cpoint> cp;
vector<Cpoint*> a;
Cpoint *root;
int n;
void init_Hull(Cpoint *a,Cpoint *b,Cpoint *c)
{
Cface *pf1=new Cface(a,b,c);
a->made.push_back(pf1);
Cface *pf2=new Cface(c,b,a);
a->made.push_back(pf2);
//pf1,pf2为相反的半三维空间(半超平面?)
pf2->set_line_back(pf1->line+1,pf1->line+0,pf1->line+2);//c,b,a
vector<Cpoint*> all;
for (int i=0;i<n;i++)
if (&cp[i]!=this->a[0]&&&cp[i]!=this->a[1]&&&cp[i]!=this->a[2])
all.push_back(&cp[i]);//非开始的确定平面的三个点
pf1->find_see(all);
pf2->find_see(all);
}
void insert_point(Cpoint *a)
{
Cline* ps=NULL;
for (auto f:a->see)//寻找加入点的可视面
if (!f->killed)
{
f->killed=a;
if (!ps)
for (int i=0;i<3;i++)
if (!see_face(f->line[i].back->f->f,a->p))
ps=f->line+i;
}
if (!ps) return;//点在凸包内
Cface* pre=NULL;
Cline* pt=ps;
while (1)
{
Cface *nf=new Cface(a,pt->back->p,pt->p);//新建一个包含新点与对应边pt的面
a->made.push_back(nf);
nf->find_see(pt->f->seen,pt->back->f->seen);//新的可见点集可能存在于pt原来的两个面的可见点集中
if (!pre)
nf->set_line_back(nf->line+2,pt->back,nf->line);
else
nf->set_line_back(pre->line+2,pt->back,pre->line[2].back);
//设置新的面相关的边所连接的两个面
pre=nf;
pt=pt->next_line();
while (pt->back->f->killed) pt=pt->back->next_line();
//寻找下一条只剩一个面的直线
if (pt==ps) break;
}
}
void clear()
{
for (int i=0;i<n;i++)
{
cp[i].made.clear();
cp[i].see.clear();
}
}
int in_Hull(point p,Cpoint *s)
{
for (auto i:s->made)
{
if (face_point_dis(i->f,p)>eps)
if (i->killed) return in_Hull(p,i->killed);//面已经被删除
else return 0;
}
return 1;
}
int in_Hull(point p)
{
return in_Hull(p,root);
}
void build(vector<point> points)
{
n=points.size();
cp.resize(n),a.resize(n);
for (int i=0;i<n;i++)
{
a[i]=&cp[i];
cp[i].p=points[i];
}
random_shuffle(a.begin(),a.end());//随机加入点
root=a[0];
init_Hull(a[0],a[1],a[2]);
for (int i=3;i<n;i++) insert_point(a[i]);
}
vector<face> get_all_face()
{
vector<face> ret;
for (auto &p:cp)
for (auto &f:p.made)
if (!f->killed) ret.push_back(f->f);//面尚未被删除
return ret;
}
point get_center()//求中心
{
point ret(0,0,0);
for (auto &p:cp) ret=ret+p.p;
return ret/cp.size();
}
}hull,Hull;
int read()
{
int x=1;
char ch;
while (ch=gc,ch<'0'||ch>'9') if (ch=='-') x=-1;
int s=ch-'0';
while (ch=gc,ch>='0'&&ch<='9') s=s*10+ch-'0';
return s*x;
}
int n,m;
int main()
{
srand(20010527);
n=read();
for (int i=1;i<=n;i++)
{
double a=read(),b=read(),c=read(),d=read();
points.push_back(point(a,b,c)/(a+b+c+d));
}
hull.build(points);
vector<face> fs=hull.get_all_face();
point center=hull.get_center();
points.clear();
for (auto &f:fs) points.push_back(face_to_point(f,center));
Hull.build(points);
m=read();
while (m--)
{
int op=read();
double a=read(),b=read(),c=read(),d=read();
int ret=0;
if (op==1) ret=hull.in_Hull(point(a,b,c)/(a+b+c+d));
else
{
face f;
f.normal=point(a-d,b-d,c-d);
f.d=-d;
if (!see_face(f,center))//如果中心满足显然可以
ret=1-Hull.in_Hull(face_to_point(f,center));
else ret=1;
}
puts(ret?"Y":"N");
}
return 0;
}
-------------本文结束感谢您的阅读-------------

本文标题:Uoj#44.【清华集训2014】Breaking Bomber

文章作者:wzf2000

发布时间:2018年04月18日 - 15:04

最后更新:2018年04月27日 - 11:04

原始链接:https://wzf2000.github.io/2018/04/18/Uoj44/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。