Skip to content

Commit 8971c03

Browse files
ChilleeChillee
committed
some debugging stuff
1 parent b6f7c82 commit 8971c03

File tree

10 files changed

+287
-43
lines changed

10 files changed

+287
-43
lines changed

‎content/geometry/HalfPlane.h‎

Lines changed: 26 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -27,18 +27,40 @@ bool cmp(Line a, Line b) {
2727
}
2828
vector<P> halfPlaneIntersection(vector<Line> vs) {
2929
sort(all(vs), cmp);
30+
// for (auto l: vs) {
31+
// cout<<l[0]<<" -> "<<l[1]<<endl;
32+
// }
3033
vector<Line> deq(sz(vs) + 5);
3134
vector<P> ans(sz(vs) + 5);
3235
deq[0] = vs[0];
3336
int ah = 0, at = 0, n = sz(vs);
3437
rep(i,1,n+1) {
3538
if (i == n) vs.push_back(deq[ah]);
36-
if (angDiff(vs[i], vs[i - 1]) == 0) continue;
37-
while (ah<at && sideOf(sp(vs[i]),ans[at-1]) < 0) at--;
38-
while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah])<0) ah++;
39+
if (angDiff(vs[i], vs[i - 1]) == 0) {
40+
continue;
41+
}
42+
43+
while (ah<at && sideOf(sp(vs[i]),ans[at-1], 1e-8) < 0) at--;
44+
while (i!=n && ah<at && sideOf(sp(vs[i]),ans[ah], 1e-8)<0) ah++;
3945
auto res = lineInter(sp(vs[i]), sp(deq[at]));
40-
if (res.first != 1) continue;
46+
if (res.first != 1) {
47+
// cout<<"46"<<endl;
48+
continue;
49+
}
4150
ans[at++] = res.second, deq[at] = vs[i];
51+
// cout<<"i: "<<i<<endl;
52+
// cout<<"ans: ";
53+
// cout<<ah<<' '<<at<<endl;
54+
// for (int j=ah; j<at; j++) {
55+
// cout<<ans[j]<<" ";
56+
// }
57+
// cout<<endl;
58+
// cout<<"deq: ";
59+
// for(int j=ah; j<=at; j++) {
60+
// cout<<deq[j][0]<<"-> "<<deq[j][1]<<' ';
61+
// }
62+
// cout<<endl;
63+
// cout<<endl;
4264
}
4365
if (at - ah <= 2) return {};
4466
return {ans.begin() + ah, ans.begin() + at};

‎content/geometry/PolygonArea.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,8 @@ template<class T>
1515
T polygonArea2(vector<Point<T>> v) {
1616
if (sz(v) < 3) return 0;
1717
T a = v.back().cross(v[0]);
18-
rep(i,0,sz(v)-1) a += v[i].cross(v[i+1]);
18+
rep(i,0,sz(v)-1) {
19+
a += v[i].cross(v[i+1]);
20+
}
1921
return a;
2022
}

‎content/geometry/lineIntersection.h‎

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -25,9 +25,18 @@ Products of three coordinates are used in intermediate steps so watch out for ov
2525

2626
template<class P>
2727
pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
28-
auto d = (e1 - s1).cross(e2 - s2);
29-
if (d == 0) // if parallel
30-
return {-(s1.cross(e1, s2) == 0), P(0, 0)};
31-
auto p = s2.cross(e1, e2), q = s2.cross(e2, s1);
32-
return {1, (s1 * p + e1 * q) / d};
28+
auto d = (e1-s1).cross(e2-s2);
29+
if (d == 0) //if parallel
30+
return {-((e1-s1).cross(s2-s1)==0 || s2==e2), P(0,0)};
31+
else
32+
return {1, s2-(e2-s2)*(e1-s1).cross(s2-s1)/d};
3333
}
34+
35+
// template<class P>
36+
// pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
37+
// auto d = (e1 - s1).cross(e2 - s2);
38+
// if (d == 0) // if parallel
39+
// return {-(s1.cross(e1, s2) == 0), P(0, 0)};
40+
// auto p = s2.cross(e1, e2), q = s2.cross(e2, s1);
41+
// return {1, (s1 * p + e1 * q) / d};
42+
// }

‎content/geometry/sideOf.h‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,5 +19,11 @@ template<class P>
1919
int sideOf(const P& s, const P& e, const P& p, double eps) {
2020
auto a = (e-s).cross(p-s);
2121
double l = (e-s).dist()*eps;
22-
return (a > l) - (a < -l);
22+
if (a<=l && a>=-l) {
23+
return 0;
24+
} else if (a > l) {
25+
return 1;
26+
} else {
27+
return -1;
28+
}
2329
}

‎content/strings/Hashing.h‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -16,10 +16,12 @@
1616
struct H {
1717
typedef uint64_t ull;
1818
ull x; H(ull x=0) : x(x) {}
19-
#define OP(O,A,B) H operator O(H o) { ull r = x; asm \
20-
(A "addq %%rdx, %0\n adcq $0,%0" : "+a"(r) : B); return r; }
21-
OP(+,,"d"(o.x)) OP(*,"mul %1\n", "r"(o.x) : "rdx")
22-
H operator-(H o) { return *this + ~o.x; }
19+
H operator+(H o) { return x + o.x + (x + o.x < x); }
20+
H operator*(H o) {
21+
ull r = x;
22+
asm("mul %1\n addq %%rdx, %0\n adcq $0,%0" : "+a"(r) : "r"(o.x) : "rdx");
23+
return r;
24+
} H operator-(H o) { return *this + ~o.x; }
2325
ull get() const { return x + !~x; }
2426
bool operator==(H o) const { return get() == o.get(); }
2527
bool operator<(H o) const { return get() < o.get(); }

‎fuzz-tests/geometry/HalfPlane‎

1000 KB
Binary file not shown.

‎fuzz-tests/geometry/HalfPlane.cpp‎

Lines changed: 128 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -302,22 +302,37 @@ void testEmpty() {
302302
assert(sz(res) == 0);
303303
}
304304
void testRandom() {
305-
int lim = 3;
306-
for (int i = 0; i < 1000000; i++) {
305+
int lim = 5;
306+
for (int i = 0; i < 10000000; i++) {
307307
vector<Line> t;
308-
for (int i = 0; i < 6; i++) {
308+
for (int i = 0; i < 4; i++) {
309309
Line cand{P(0, 0), P(0, 0)};
310310
while (cand[0] == cand[1])
311311
cand = {randPt(lim), randPt(lim)};
312+
Point<int> s1 = Point<int>(cand[0].x, cand[0].y);
313+
Point<int> e1 = Point<int>(cand[1].x, cand[1].y);
314+
bool fail = false;
315+
for (auto j: t) {
316+
Point<int> s2 = Point<int>(j[0].x, j[0].y);
317+
Point<int> e2 = Point<int>(j[1].x, j[1].y);
318+
if (lineInter(s1, e1, s2, e2).first == -1) {
319+
fail = false;
320+
break;
321+
}
322+
}
323+
if (fail) {
324+
i--;
325+
continue;
326+
}
312327
t.push_back(cand);
313328
}
314329
addInf(t, lim);
315330
auto res = halfPlaneIntersection(t);
316331
double area = sjtu::halfPlaneIntersection(t);
317-
// double resArea = polygonArea2(res)/2;
318-
double resArea = mit::Intersection_Area(convert(t));
332+
double resArea = abs(polygonArea2(res) / 2);
333+
// double resArea = mit::Intersection_Area(convert(t));
319334
double diff = abs(area - resArea);
320-
if (diff > EPS) {
335+
if (diff > EPS || isnan(diff)) {
321336
cout << diff << ' ' << area << ' ' << resArea << endl;
322337
for (auto i : t)
323338
cout << i[0] << "->" << i[1] << ' ';
@@ -330,36 +345,121 @@ void testRandom() {
330345
cout << i << ',';
331346
cout << endl;
332347
}
333-
assert(diff < EPS);
348+
assert(diff <= EPS);
334349
}
335350
}
351+
vector<P> convexHull(vector<P> pts) {
352+
if (sz(pts) <= 1) return pts;
353+
sort(all(pts));
354+
vector<P> h(sz(pts)+1);
355+
int s = 0, t = 0;
356+
for (int it = 2; it--; s = --t, reverse(all(pts)))
357+
trav(p, pts) {
358+
while (t >= s + 2 && h[t-2].cross(h[t-1], p) <= 0) t--;
359+
h[t++] = p;
360+
}
361+
return {h.begin(), h.begin() + t - (t == 2 && h[0] == h[1])};
362+
}
363+
336364

337365
int main() {
338-
test1();
339-
testInf();
340-
testLine();
341-
testPoint();
342-
testEmpty();
343-
testRandom();
366+
// test1();
367+
// testInf();
368+
// testLine();
369+
// testPoint();
370+
// testRandom();
371+
// return 0;
372+
// testEmpty();
344373
// Case that messes with precision
345-
vector<Line> t({{P(8, 9), P(8, 2)},
346-
{P(3, 9), P(5, 2)},
347-
{P(8, 2), P(8, 3)},
348-
{P(7, 2), P(1, 7)},
349-
{P(1, 0), P(7, 1)},
350-
{P(9, 2), P(5, 6)},
351-
{P(10, 10), P(-10, 10)},
352-
{P(-10, 10), P(-10, -10)},
353-
{P(-10, -10), P(10, -10)},
354-
{P(10, -10), P(10, 10)}});
355-
auto res = halfPlaneIntersection(t);
356-
cout << polygonArea2(res) << endl;
357-
cout << sjtu::halfPlaneIntersection(t) << endl;
358-
cout << mit::Intersection_Area(convert(t))<<endl;
374+
// vector<Line> cases({{P(1,0),P(0,2)},{P(0,1),P(2,0)},{P(0,0),P(1,1)},{P(2,2),P(1,1)},{P(3,3),P(-3,3)},{P(-3,3),P(-3,-3)},{P(-3,-3),P(3,-3)},{P(3,-3),P(3,3)}});
375+
376+
// auto pts = halfPlaneIntersection(cases);
377+
// for (auto p: pts) {
378+
// cout<<p<<' ';
379+
// }
380+
// cout<<endl;
381+
// cout<<polygonArea2(pts)<<endl;
382+
// return 0;
383+
384+
array<int, 3> mn = {1e6, 1e6, 10};
385+
for (int i = 0; i < 10000000; i++) {
386+
ll offset = rand() % ll(1e9);
387+
ll mul = rand() % ll(1e9);
388+
vector<pair<int, vector<Line>>> cases;
389+
390+
cases.push_back({10,
391+
{{P(offset + mul * 8, offset + mul * 9), P(offset + mul * 8, offset + mul * 2)},
392+
{P(offset + mul * 3, offset + mul * 9), P(offset + mul * 5, offset + mul * 2)},
393+
{P(offset + mul * 8, offset + mul * 2), P(offset + mul * 8, offset + mul * 3)},
394+
{P(offset + mul * 7, offset + mul * 2), P(offset + mul * 1, offset + mul * 7)},
395+
{P(offset + mul * 1, offset + mul * 0), P(offset + mul * 7, offset + mul * 1)},
396+
{P(offset + mul * 9, offset + mul * 2), P(offset + mul * 5, offset + mul * 6)},
397+
{P(offset + mul * 10, offset + mul * 10), P(offset + mul * -10, offset + mul * 10)},
398+
{P(offset + mul * -10, offset + mul * 10), P(offset + mul * -10, offset + mul * -10)},
399+
{P(offset + mul * -10, offset + mul * -10), P(offset + mul * 10, offset + mul * -10)},
400+
{P(offset + mul * 10, offset + mul * -10), P(offset + mul * 10, offset + mul * 10)}}});
401+
cases.push_back({5,
402+
{{P(offset + mul * 0, offset + mul * 1), P(offset + mul * 4, offset + mul * 0)},
403+
{P(offset + mul * 0, offset + mul * 2), P(offset + mul * 2, offset + mul * 0)},
404+
{P(offset + mul * 1, offset + mul * 0), P(offset + mul * 3, offset + mul * 4)},
405+
{P(offset + mul * 4, offset + mul * 2), P(offset + mul * 0, offset + mul * 0)},
406+
{P(offset + mul * 4, offset + mul * 2), P(offset + mul * 2, offset + mul * 2)},
407+
{P(offset + mul * 0, offset + mul * 3), P(offset + mul * 1, offset + mul * 1)},
408+
{P(offset + mul * 5, offset + mul * 5), P(offset + mul * -5, offset + mul * 5)},
409+
{P(offset + mul * -5, offset + mul * 5), P(offset + mul * -5, offset + mul * -5)},
410+
{P(offset + mul * -5, offset + mul * -5), P(offset + mul * 5, offset + mul * -5)},
411+
{P(offset + mul * 5, offset + mul * -5), P(offset + mul * 5, offset + mul * 5)}}});
412+
cases.push_back({20,
413+
{{P(offset + mul * 10, offset + mul * 12), P(offset + mul * 16, offset + mul * 5)},
414+
{P(offset + mul * 10, offset + mul * 12), P(offset + mul * 4, offset + mul * 19)},
415+
{P(offset + mul * 18, offset + mul * 15), P(offset + mul * 17, offset + mul * 7)},
416+
{P(offset + mul * 12, offset + mul * 13), P(offset + mul * 14, offset + mul * 6)},
417+
{P(offset + mul * 16, offset + mul * 3), P(offset + mul * 4, offset + mul * 11)},
418+
{P(offset + mul * 12, offset + mul * 8), P(offset + mul * 9, offset + mul * 0)},
419+
{P(offset + mul * 20, offset + mul * 20), P(offset + mul * -20, offset + mul * 20)},
420+
{P(offset + mul * -20, offset + mul * 20), P(offset + mul * -20, offset + mul * -20)},
421+
{P(offset + mul * -20, offset + mul * -20), P(offset + mul * 20, offset + mul * -20)},
422+
{P(offset + mul * 20, offset + mul * -20), P(offset + mul * 20, offset + mul * 20)}}});
423+
// {P(5,8),P(4,9)},{P(9,1),P(5,8)},{P(9,7),P(7,6)},{P(3,7),P(8,5)},{P(6,6),P(8,4)},{P(6,1),P(7,2)},{P(10,10),P(-10,10)},{P(-10,10),P(-10,-10)},{P(-10,-10),P(10,-10)},{P(10,-10),P(10,10)},
424+
int idx = 0;
425+
for (auto tmp : cases) {
426+
auto t = tmp.second;
427+
auto lim = tmp.first;
428+
auto res = halfPlaneIntersection(t);
429+
auto ours = polygonArea2(res);
430+
if (abs(ours) > 1e-3) {
431+
if (mn[0] + mn[1] * mn[2] >= offset + mul * tmp.first) {
432+
cout << "case: " << idx << endl;
433+
cout << offset << ' ' << mul << ' ' << offset + mul * tmp.first << endl;
434+
cout << ours << ' ' << sjtu::halfPlaneIntersection(t) << endl;
435+
436+
double mx = 0;
437+
for (auto i: res) {
438+
for (auto j: res)
439+
mx = max(mx, (i - j).dist());
440+
}
441+
cout<<"mx: "<<mx<<endl;
442+
for (auto i: res) {
443+
cout<<i<<' ';
444+
}
445+
cout<<endl;
446+
cout << endl;
447+
mn = {offset, mul, tmp.first};
448+
}
449+
// cout << ours << endl;
450+
// cout << sjtu::halfPlaneIntersection(t) << endl;
451+
// auto mits = mit::Intersection_Area(convert(t));
452+
// cout << mits<<endl;
453+
// cout << endl;
454+
}
455+
idx++;
456+
}
457+
}
458+
cout << "mn: " << mn[0] << ' ' << mn[1] << endl;
459+
cout << mn[0] + mn[1] * mn[2] << endl;
359460
// Failure case for mit's half plane cod
360461
// vector<Line> t({{P(9, 8), P(9, 1)}, {P(3, 3), P(3, 5)}, {P(7, 6), P(0, 8)}});
361462
// Failure case for old code.
362-
// vector<Line> t({{P(3,0),P(3,3)},{P(5,3), P(5,0)}, {P(4,-2), P(0,1)}, {P(3,-1), P(0,-1)}});
363463
// addInf(t);
364464
// cout << polygonArea2(halfPlaneIntersection(t)) << endl;
365465
// addInf(t);

‎fuzz-tests/geometry/t‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
(2,2) -> (1,1)
2+
(-3,3) -> (-3,-3)
3+
(0,1) -> (2,0)
4+
(-3,-3) -> (3,-3)
5+
(0,0) -> (1,1)
6+
(3,-3) -> (3,3)
7+
(1,0) -> (0,2)
8+
(3,3) -> (-3,3)
9+
i: 1
10+
ans: 0 1
11+
(-3,-3)
12+
deq: (2,2)-> (1,1) (-3,3)-> (-3,-3)
13+
14+
i: 2
15+
ans: 0 1
16+
(0.666667,0.666667)
17+
deq: (2,2)-> (1,1) (0,1)-> (2,0)
18+
19+
i: 3
20+
ans: 0 2
21+
(0.666667,0.666667) (8,-3)
22+
deq: (2,2)-> (1,1) (0,1)-> (2,0) (-3,-3)-> (3,-3)
23+
24+
i: 4
25+
ans: 0 2
26+
(0.666667,0.666667) (0.666667,0.666667)
27+
deq: (2,2)-> (1,1) (0,1)-> (2,0) (0,0)-> (1,1)
28+
29+
i: 5
30+
ans: 0 3
31+
(0.666667,0.666667) (0.666667,0.666667) (3,3)
32+
deq: (2,2)-> (1,1) (0,1)-> (2,0) (0,0)-> (1,1) (3,-3)-> (3,3)
33+
34+
i: 6
35+
ans: 1 3
36+
(0.666667,0.666667) (0.666667,0.666667)
37+
deq: (0,1)-> (2,0) (0,0)-> (1,1) (1,0)-> (0,2)
38+
39+
i: 7
40+
ans: 1 4
41+
(0.666667,0.666667) (0.666667,0.666667) (-0.5,3)
42+
deq: (0,1)-> (2,0) (0,0)-> (1,1) (1,0)-> (0,2) (3,3)-> (-3,3)
43+
44+
i: 8
45+
ans: 1 5
46+
(0.666667,0.666667) (0.666667,0.666667) (-0.5,3) (-4,3)
47+
deq: (0,1)-> (2,0) (0,0)-> (1,1) (1,0)-> (0,2) (3,3)-> (-3,3) (0,1)-> (2,0)
48+
49+
(0.666667,0.666667) (0.666667,0.666667) (-0.5,3) (-4,3)
50+
8.16667

0 commit comments

Comments
 (0)