diff --git a/stress-tests/data-structures/LineContainer.cpp b/stress-tests/data-structures/LineContainer.cpp index 1199349a6..d810cca5c 100644 --- a/stress-tests/data-structures/LineContainer.cpp +++ b/stress-tests/data-structures/LineContainer.cpp @@ -6,38 +6,38 @@ namespace other { // source: https://github.com/niklasb/contest-algos/blob/master/convex_hull/dynamic.cpp const ll is_query = -(1LL<<62); struct Line { - ll m, b; - mutable function succ; - bool operator<(const Line& rhs) const { - if (rhs.b != is_query) return m < rhs.m; - const Line* s = succ(); - if (!s) return 0; - ll x = rhs.m; - return b - s->b < (s->m - m) * x; - } + ll m, b; + mutable function succ; + bool operator<(const Line& rhs) const { + if (rhs.b != is_query) return m < rhs.m; + const Line* s = succ(); + if (!s) return 0; + ll x = rhs.m; + return b - s->b < (s->m - m) * x; + } }; struct HullDynamic : public multiset { // will maintain upper hull for maximum - bool bad(iterator y) { - auto z = next(y); - if (y == begin()) { - if (z == end()) return 0; - return y->m == z->m && y->b <= z->b; - } - auto x = prev(y); - if (z == end()) return y->m == x->m && y->b <= x->b; - return (x->b - y->b)*(z->m - y->m) >= (y->b - z->b)*(y->m - x->m); - } - void add(ll m, ll b) { - auto y = insert({ m, b }); - y->succ = [=] { return next(y) == end() ? 0 : &*next(y); }; - if (bad(y)) { erase(y); return; } - while (next(y) != end() && bad(next(y))) erase(next(y)); - while (y != begin() && bad(prev(y))) erase(prev(y)); - } - ll query(ll x) { - auto l = *lower_bound((Line) { x, is_query }); - return l.m * x + l.b; - } + bool bad(iterator y) { + auto z = next(y); + if (y == begin()) { + if (z == end()) return 0; + return y->m == z->m && y->b <= z->b; + } + auto x = prev(y); + if (z == end()) return y->m == x->m && y->b <= x->b; + return (x->b - y->b)*(z->m - y->m) >= (y->b - z->b)*(y->m - x->m); + } + void add(ll m, ll b) { + auto y = insert({ m, b }); + y->succ = [=] { return next(y) == end() ? 0 : &*next(y); }; + if (bad(y)) { erase(y); return; } + while (next(y) != end() && bad(next(y))) erase(next(y)); + while (y != begin() && bad(prev(y))) erase(prev(y)); + } + ll query(ll x) { + auto l = *lower_bound((Line) { x, is_query }); + return l.m * x + l.b; + } }; } diff --git a/stress-tests/geometry/CircleLine.cpp b/stress-tests/geometry/CircleLine.cpp index ee322252a..cfe9d8609 100644 --- a/stress-tests/geometry/CircleLine.cpp +++ b/stress-tests/geometry/CircleLine.cpp @@ -6,20 +6,20 @@ typedef Point P; int main() { - { - auto res = circleLine(P(0, 0), 1, P(-1, -1), P(1, 1)); - assert(res.size() == 2); - assert((res[1]-P(sqrt(2)/2, sqrt(2)/2)).dist() < 1e-8); - } - { - auto res = circleLine(P(0, 0), 1, P(-5, 1), P(5, 1)); - assert(res.size() == 1); - assert((res[0]-P(0,1)).dist() < 1e-8); - } - { - auto res = circleLine(P(4, 4), 1, P(0, 0), P(5, 0)); - assert(res.size() == 0); - } + { + auto res = circleLine(P(0, 0), 1, P(-1, -1), P(1, 1)); + assert(res.size() == 2); + assert((res[1]-P(sqrt(2)/2, sqrt(2)/2)).dist() < 1e-8); + } + { + auto res = circleLine(P(0, 0), 1, P(-5, 1), P(5, 1)); + assert(res.size() == 1); + assert((res[0]-P(0,1)).dist() < 1e-8); + } + { + auto res = circleLine(P(4, 4), 1, P(0, 0), P(5, 0)); + assert(res.size() == 0); + } rep(it,0,100000) { P a = randIntPt(5); P b = randIntPt(5); @@ -46,5 +46,5 @@ int main() { assert(!points.empty()); } } - cout<<"Tests passed!"< P; long double areaCT(P pa, P pb, long double r) { - if (pa.dist() < pb.dist()) swap(pa, pb); - if (sgn(pb.dist()) == 0) return 0; - long double a = pb.dist(), b = pa.dist(), c = (pb - pa).dist(); - long double sinB = fabs(pb.cross(pb - pa) / a / c), cosB = pb.dot(pb - pa) / a / c, - sinC = fabs(pa.cross(pb) / a / b), cosC = pa.dot(pb) / a / b; - long double B = atan2(sinB, cosB), C = atan2(sinC, cosC); - if (a > r) { - long double S = C / 2 * r * r, h = a * b * sinC / c; - if (h < r && B < M_PI / 2) - S -= (acos(h / r) * r * r - h * sqrt(r * r - h * h)); - return S; - } else if (b > r) { - long double theta = M_PI - B - asin(sinB / r * a); - return a * r * sin(theta) / 2 + (C - theta) / 2 * r * r; - } else return sinC * a * b / 2; + if (pa.dist() < pb.dist()) swap(pa, pb); + if (sgn(pb.dist()) == 0) return 0; + long double a = pb.dist(), b = pa.dist(), c = (pb - pa).dist(); + long double sinB = fabs(pb.cross(pb - pa) / a / c), cosB = pb.dot(pb - pa) / a / c, + sinC = fabs(pa.cross(pb) / a / b), cosC = pa.dot(pb) / a / b; + long double B = atan2(sinB, cosB), C = atan2(sinC, cosC); + if (a > r) { + long double S = C / 2 * r * r, h = a * b * sinC / c; + if (h < r && B < M_PI / 2) + S -= (acos(h / r) * r * r - h * sqrt(r * r - h * h)); + return S; + } else if (b > r) { + long double theta = M_PI - B - asin(sinB / r * a); + return a * r * sin(theta) / 2 + (C - theta) / 2 * r * r; + } else return sinC * a * b / 2; } long double circlePoly(P c, long double r, vector

poly) { - long double area = 0; - rep(i,0,sz(poly)){ - auto a = poly[i] - c, b = poly[(i+1)%sz(poly)] - c; - area += areaCT(a, b, r) * sgn(a.cross(b)); - } - return area; + long double area = 0; + rep(i,0,sz(poly)){ + auto a = poly[i] - c, b = poly[(i+1)%sz(poly)] - c; + area += areaCT(a, b, r) * sgn(a.cross(b)); + } + return area; } } signed main() { - ios::sync_with_stdio(0); - cin.tie(0); - const int lim=5; - for (int i=0;i<100000; i++) { + ios::sync_with_stdio(0); + cin.tie(0); + const int lim=5; + for (int i=0;i<100000; i++) { - vector> pts; - for (int j=0; j<10; j++) { - int x = rand()%lim, y = rand()%lim; - pts.push_back(Point(x, y)); - } + vector> pts; + for (int j=0; j<10; j++) { + int x = rand()%lim, y = rand()%lim; + pts.push_back(Point(x, y)); + } - auto polyInt = genPolygon(pts); + auto polyInt = genPolygon(pts); - int cx = rand()%lim, cy = rand()%lim; - auto c = P(cx, cy); - auto c2 = orig::P(cx, cy); - double r= rand()%(2*lim); + int cx = rand()%lim, cy = rand()%lim; + auto c = P(cx, cy); + auto c2 = orig::P(cx, cy); + double r= rand()%(2*lim); - vector

poly; - vector poly2; - for (auto j: polyInt) { - poly.push_back(P(j.x, j.y)); - poly2.push_back(orig::P(j.x, j.y)); - } - auto res1 = circlePoly(c, r, poly); - auto res2 = orig::circlePoly(c2, r, poly2); + vector

poly; + vector poly2; + for (auto j: polyInt) { + poly.push_back(P(j.x, j.y)); + poly2.push_back(orig::P(j.x, j.y)); + } + auto res1 = circlePoly(c, r, poly); + auto res2 = orig::circlePoly(c2, r, poly2); - if (abs(res1 - res2) > 1e-8) { - cout< 1e-8) { + cout< P; signed main() { - for (int i = 0; i < 1000000; i++) { - P c1 = randIntPt(5), c2 = randIntPt(5); - double r1 = sqrt(rand()%20), r2 = sqrt(rand()%20); - for (auto sgn : {-1, 1}) { - auto tans = tangents(c1, r1, c2, sgn * r2); + for (int i = 0; i < 1000000; i++) { + P c1 = randIntPt(5), c2 = randIntPt(5); + double r1 = sqrt(rand()%20), r2 = sqrt(rand()%20); + for (auto sgn : {-1, 1}) { + auto tans = tangents(c1, r1, c2, sgn * r2); - if (tans.size() ==1) { - assert((tans[0].first - tans[0].second).dist() < 1e-8); - assert(abs((tans[0].first-c1).dist() - r1) < 1e-8); - assert(abs((tans[0].first-c2).dist() - r2) < 1e-8); - } else if (tans.size() == 2) { - for (auto l : tans) { - assert(abs(abs(lineDist(l.first, l.second, c1))-r1) < 1e-8); - assert(abs(abs(lineDist(l.first, l.second, c2))-r2) < 1e-8); - } - } - } - } - cout<<"Tests passed!"<& S) { } int main() { - const int SZ = 1e2; - rep(t,0,100000) { - const int GRID=1e3; - vector

pts(SZ); - rep(i,0,SZ) pts[i] = P(rand()%GRID, rand()%GRID); - auto res = convexHull(pts); - auto res2 = old::convexHull(pts); - assert(sz(res) == sz(res2)); - rep(i,0,sz(res2)) { - assert(pts[res2[i]] == res[i]); - } - } - cout<<"Tests passed!"< pts(SZ); + rep(i,0,SZ) pts[i] = P(rand()%GRID, rand()%GRID); + auto res = convexHull(pts); + auto res2 = old::convexHull(pts); + assert(sz(res) == sz(res2)); + rep(i,0,sz(res2)) { + assert(pts[res2[i]] == res[i]); + } + } + cout<<"Tests passed!"< P; int main() { - cin.sync_with_stdio(0); - cin.tie(0); - const int lim = 5; - for (int i = 0; i < 100000; i++) { - P p = P(rand() % lim, rand() % lim); - P a = P(rand() % lim, rand() % lim); - P b = P(rand() % lim, rand() % lim); - while (a == b) - b = P(rand() % lim, rand() % lim); - auto proj = lineProj(a, b, p, false); - auto refl = lineProj(a, b, p, true); - assert(lineDist(a, b, proj) < 1e-8); - auto manProj = (refl + p) / 2; - assert((proj-manProj).dist() < 1e-8); - } - cout<<"Tests passed!"< P; typedef int T; T rectilinear_mst_n(vector

ps) { - struct edge { int src, dst; T weight; }; - vector edges; + struct edge { int src, dst; T weight; }; + vector edges; - auto dist = [&](int i, int j) { - return abs((ps[i]-ps[j]).x) + abs((ps[i]-ps[j]).y); - }; - for (int i = 0; i < sz(ps); ++i) - for (int j = i+1; j < sz(ps); ++j) - edges.push_back({i, j, dist(i,j)}); - T cost = 0; - sort(all(edges), [](edge a, edge b) { return a.weight < b.weight; }); - UF uf(sz(ps)); - for (auto e: edges) - if (uf.join(e.src, e.dst)) - cost += e.weight; - return cost; + auto dist = [&](int i, int j) { + return abs((ps[i]-ps[j]).x) + abs((ps[i]-ps[j]).y); + }; + for (int i = 0; i < sz(ps); ++i) + for (int j = i+1; j < sz(ps); ++j) + edges.push_back({i, j, dist(i,j)}); + T cost = 0; + sort(all(edges), [](edge a, edge b) { return a.weight < b.weight; }); + UF uf(sz(ps)); + for (auto e: edges) + if (uf.join(e.src, e.dst)) + cost += e.weight; + return cost; } signed main() { - for (int t=0; t<10000; t++) { - const int max_coord = rand() % 300 + 1; - const int num_pts = rand() % 100; - vector

pts; - for (int i = 0; i < num_pts; ++i) { - int x = rand() % max_coord - max_coord / 2; - int y = rand() % max_coord - max_coord / 2; - pts.push_back(P(x,y)); - } - auto edges = manhattanMST(pts); - assert(edges.size() <= 4*pts.size()); - sort(all(edges)); - UF uf(sz(pts)); - int cost = 0, joined = 0; - for (auto e: edges) if (uf.join(e[1], e[2])) cost += e[0], joined++; - if (num_pts > 0) assert(joined == num_pts - 1); - assert(cost == rectilinear_mst_n(pts)); - } - cout<<"Tests passed!"< pts; + for (int i = 0; i < num_pts; ++i) { + int x = rand() % max_coord - max_coord / 2; + int y = rand() % max_coord - max_coord / 2; + pts.push_back(P(x,y)); + } + auto edges = manhattanMST(pts); + assert(edges.size() <= 4*pts.size()); + sort(all(edges)); + UF uf(sz(pts)); + int cost = 0, joined = 0; + for (auto e: edges) if (uf.join(e[1], e[2])) cost += e[0], joined++; + if (num_pts > 0) assert(joined == num_pts - 1); + assert(cost == rectilinear_mst_n(pts)); + } + cout<<"Tests passed!"< eps) - (x < -eps); } @@ -37,67 +37,67 @@ db vect(pt p1, pt p2) { return p1.x * p2.y - p1.y * p2.x; } db scal(pt p1, pt p2) { return p1.x * p2.x + p1.y * p2.y; } db polygon_union(vector poly[], int n) { - auto ratio = [](pt A, pt B, pt O) { - return !sgn(A.x - B.x) ? (O.y - A.y) / (B.y - A.y) : (O.x - A.x) / (B.x - A.x); - }; - db ret = 0; - for (int i = 0; i < n; ++i) { - for (size_t v = 0; v < poly[i].size(); ++v) { - pt A = poly[i][v], B = poly[i][(v + 1) % poly[i].size()]; - vector> segs; - segs.emplace_back(0, 0), segs.emplace_back(1, 0); - for (int j = 0; j < n; ++j) - if (i != j) { - for (size_t u = 0; u < poly[j].size(); ++u) { - pt C = poly[j][u], D = poly[j][(u + 1) % poly[j].size()]; - int sc = sgn(vect(B - A, C - A)), sd = sgn(vect(B - A, D - A)); - if (!sc && !sd) { - if (sgn(scal(B - A, D - C)) > 0 && i > j) { - segs.emplace_back(ratio(A, B, C), 1), segs.emplace_back(ratio(A, B, D), -1); - } - } else { - db sa = vect(D - C, A - C), sb = vect(D - C, B - C); - if (sc >= 0 && sd < 0) - segs.emplace_back(sa / (sa - sb), 1); - else if (sc < 0 && sd >= 0) - segs.emplace_back(sa / (sa - sb), -1); - } - } - } - sort(segs.begin(), segs.end()); - db pre = min(max(segs[0].first, 0.0), 1.0), now, sum = 0; - int cnt = segs[0].second; - for (size_t j = 1; j < segs.size(); ++j) { - now = min(max(segs[j].first, 0.0), 1.0); - if (!cnt) - sum += now - pre; - cnt += segs[j].second; - pre = now; - } - ret += vect(A, B) * sum; - } - } - return ret / 2; + auto ratio = [](pt A, pt B, pt O) { + return !sgn(A.x - B.x) ? (O.y - A.y) / (B.y - A.y) : (O.x - A.x) / (B.x - A.x); + }; + db ret = 0; + for (int i = 0; i < n; ++i) { + for (size_t v = 0; v < poly[i].size(); ++v) { + pt A = poly[i][v], B = poly[i][(v + 1) % poly[i].size()]; + vector> segs; + segs.emplace_back(0, 0), segs.emplace_back(1, 0); + for (int j = 0; j < n; ++j) + if (i != j) { + for (size_t u = 0; u < poly[j].size(); ++u) { + pt C = poly[j][u], D = poly[j][(u + 1) % poly[j].size()]; + int sc = sgn(vect(B - A, C - A)), sd = sgn(vect(B - A, D - A)); + if (!sc && !sd) { + if (sgn(scal(B - A, D - C)) > 0 && i > j) { + segs.emplace_back(ratio(A, B, C), 1), segs.emplace_back(ratio(A, B, D), -1); + } + } else { + db sa = vect(D - C, A - C), sb = vect(D - C, B - C); + if (sc >= 0 && sd < 0) + segs.emplace_back(sa / (sa - sb), 1); + else if (sc < 0 && sd >= 0) + segs.emplace_back(sa / (sa - sb), -1); + } + } + } + sort(segs.begin(), segs.end()); + db pre = min(max(segs[0].first, 0.0), 1.0), now, sum = 0; + int cnt = segs[0].second; + for (size_t j = 1; j < segs.size(); ++j) { + now = min(max(segs[j].first, 0.0), 1.0); + if (!cnt) + sum += now - pre; + cnt += segs[j].second; + pre = now; + } + ret += vect(A, B) * sum; + } + } + return ret / 2; } } // namespace blackhorse namespace approximate { #include "../../content/geometry/InsidePolygon.h" double polygonUnion(vector> &polygons, int lim) { - int cnt = 0; - int total = 0; - for (double y = -lim + 1e-5; y < lim; y += lim / 500.0) { - for (double x = -lim + 1.1e-5; x < lim; x += lim / 500.0) { - total++; - for (auto &i : polygons) { - if (inPolygon(i, P(x, y))) { - cnt++; - break; - } - } - } - } - return lim * lim * 4 * cnt / double(total); + int cnt = 0; + int total = 0; + for (double y = -lim + 1e-5; y < lim; y += lim / 500.0) { + for (double x = -lim + 1.1e-5; x < lim; x += lim / 500.0) { + total++; + for (auto &i : polygons) { + if (inPolygon(i, P(x, y))) { + cnt++; + break; + } + } + } + } + return lim * lim * 4 * cnt / double(total); } } // namespace approximate @@ -113,52 +113,52 @@ inline int sgn(db x) { return (x > 1e-8) - (x < -1e-8); } typedef complex cpoi; db polygon_union(vector py[], int n) { - auto ratio = [](cpoi &a, cpoi &b, cpoi &c) { - cpoi x = b - a, y = c - a; - if (sgn(re(x)) == 0) - return im(y) / im(x); - return re(y) / re(x); - }; - db ret = 0; - for (int i = 0; i < n; ++i) - for (size_t v = 0; v < py[i].size(); ++v) { - cpoi a = py[i][v], b = py[i][(v + 1) % py[i].size()]; - vector> segs = {{0, 0}, {1, 0}}; - for (int j = 0; j < n; ++j) - if (i != j) - for (size_t u = 0; u < py[j].size(); ++u) { - cpoi c = py[j][u], d = py[j][(u + 1) % py[j].size()]; - int sc = sgn(im(conj(b - a) * (c - a))); - int sd = sgn(im(conj(b - a) * (d - a))); - if (!sc && !sd) { - if (sgn(re(conj(b - a) * (d - c))) > 0 && i > j) { - segs.pb({ratio(a, b, c), +1}); - segs.pb({ratio(a, b, d), -1}); - } - } else { - db sa = im(conj(d - c) * (a - c)); - db sb = im(conj(d - c) * (b - c)); - if (sc >= 0 && sd < 0) - segs.pb({sa / (sa - sb), 1}); - else if (sc < 0 && sd >= 0) - segs.pb({sa / (sa - sb), -1}); - } - } - sort(segs.begin(), segs.end()); - db pre = min(max(segs[0].fir, 0.0), 1.0); - db cur, sum = 0; - int cnt = segs[0].sec; - for (size_t j = 1; j < segs.size(); ++j) { - cur = min(max(segs[j].fir, 0.0), 1.0); - if (!cnt) - sum += cur - pre; - cnt += segs[j].sec; - pre = cur; - } - ret += im(conj(a) * b) * sum; - } - ret = abs(ret) * 0.5; - return ret; + auto ratio = [](cpoi &a, cpoi &b, cpoi &c) { + cpoi x = b - a, y = c - a; + if (sgn(re(x)) == 0) + return im(y) / im(x); + return re(y) / re(x); + }; + db ret = 0; + for (int i = 0; i < n; ++i) + for (size_t v = 0; v < py[i].size(); ++v) { + cpoi a = py[i][v], b = py[i][(v + 1) % py[i].size()]; + vector> segs = {{0, 0}, {1, 0}}; + for (int j = 0; j < n; ++j) + if (i != j) + for (size_t u = 0; u < py[j].size(); ++u) { + cpoi c = py[j][u], d = py[j][(u + 1) % py[j].size()]; + int sc = sgn(im(conj(b - a) * (c - a))); + int sd = sgn(im(conj(b - a) * (d - a))); + if (!sc && !sd) { + if (sgn(re(conj(b - a) * (d - c))) > 0 && i > j) { + segs.pb({ratio(a, b, c), +1}); + segs.pb({ratio(a, b, d), -1}); + } + } else { + db sa = im(conj(d - c) * (a - c)); + db sb = im(conj(d - c) * (b - c)); + if (sc >= 0 && sd < 0) + segs.pb({sa / (sa - sb), 1}); + else if (sc < 0 && sd >= 0) + segs.pb({sa / (sa - sb), -1}); + } + } + sort(segs.begin(), segs.end()); + db pre = min(max(segs[0].fir, 0.0), 1.0); + db cur, sum = 0; + int cnt = segs[0].sec; + for (size_t j = 1; j < segs.size(); ++j) { + cur = min(max(segs[j].fir, 0.0), 1.0); + if (!cnt) + sum += cur - pre; + cnt += segs[j].sec; + pre = cur; + } + ret += im(conj(a) * b) * sum; + } + ret = abs(ret) * 0.5; + return ret; } } // namespace lovelive @@ -169,59 +169,59 @@ P rndUlp(int lim, long long ulps = 5) { return P(randNearIntUlps(lim, ulps), ran P rndEps(int lim, double eps) { return P(randNearIntEps(lim, eps), randNearIntEps(lim, eps)); } void testRandom(int n, int numPts = 10, int lim = 5, bool brute = false) { - vector> polygons; - for (int i = 0; i < n; i++) { - vector

pts; - int k = randIncl(3, numPts); - for (int j = 0; j < k; j++) { - pts.push_back(randPt(lim)); // rndEps(lim, 1e-10)); - } - polygons.push_back(genPolygon(pts)); - if (polygonArea2(polygons.back()) < 0) { - reverse(all(polygons.back())); - } - } - auto val1 = polyUnion(polygons); - vector> polygons2; - for (auto i : polygons) { - vector t; - for (auto j : i) - t.push_back({j.x, j.y}); - polygons2.push_back(t); - } - vector> polygons3; - for (auto i : polygons) { - vector t; - for (auto j : i) - t.push_back({j.x, j.y}); - polygons3.push_back(t); - } - auto val3 = blackhorse::polygon_union(polygons2.data(), sz(polygons2)); - auto val4 = lovelive::polygon_union(polygons3.data(), sz(polygons3)); - if (abs(val1 - val3) > 1e-8 || abs(val1 - val4) > 1e-8) { - rep(i, 0, n) { - for (auto &x : polygons[i]) { - cout << x << ' '; - } - cout << endl; - } - abort(); - } + vector> polygons; + for (int i = 0; i < n; i++) { + vector

pts; + int k = randIncl(3, numPts); + for (int j = 0; j < k; j++) { + pts.push_back(randPt(lim)); // rndEps(lim, 1e-10)); + } + polygons.push_back(genPolygon(pts)); + if (polygonArea2(polygons.back()) < 0) { + reverse(all(polygons.back())); + } + } + auto val1 = polyUnion(polygons); + vector> polygons2; + for (auto i : polygons) { + vector t; + for (auto j : i) + t.push_back({j.x, j.y}); + polygons2.push_back(t); + } + vector> polygons3; + for (auto i : polygons) { + vector t; + for (auto j : i) + t.push_back({j.x, j.y}); + polygons3.push_back(t); + } + auto val3 = blackhorse::polygon_union(polygons2.data(), sz(polygons2)); + auto val4 = lovelive::polygon_union(polygons3.data(), sz(polygons3)); + if (abs(val1 - val3) > 1e-8 || abs(val1 - val4) > 1e-8) { + rep(i, 0, n) { + for (auto &x : polygons[i]) { + cout << x << ' '; + } + cout << endl; + } + abort(); + } } int main() { - // int s = (int)time(0); - int s = 1; - // cout << "seed " << s << endl; - srand(s); - for (int i = 0; i < 100; i++) { - testRandom(2, 5, 5); - } - for (int i = 0; i < 100; i++) { - testRandom(2, 10, 2); - } - for (int i = 0; i < 50; i++) { - testRandom(5, 100, 5); - } - cout << "Tests passed!" << endl; + // int s = (int)time(0); + int s = 1; + // cout << "seed " << s << endl; + srand(s); + for (int i = 0; i < 100; i++) { + testRandom(2, 5, 5); + } + for (int i = 0; i < 100; i++) { + testRandom(2, 10, 2); + } + for (int i = 0; i < 50; i++) { + testRandom(5, 100, 5); + } + cout << "Tests passed!" << endl; } diff --git a/stress-tests/geometry/SegmentIntersection.cpp b/stress-tests/geometry/SegmentIntersection.cpp index b633bd87c..e2e65b504 100644 --- a/stress-tests/geometry/SegmentIntersection.cpp +++ b/stress-tests/geometry/SegmentIntersection.cpp @@ -5,56 +5,56 @@ namespace oldImpl { template int segmentIntersection(const P& s1, const P& e1, - const P& s2, const P& e2, P& r1, P& r2) { - if (e1==s1) { - if (e2==s2) { - if (e1==e2) { r1 = e1; return 1; } //all equal - else return 0; //different point segments - } else return segmentIntersection(s2,e2,s1,e1,r1,r2);//swap - } - //segment directions and separation - P v1 = e1-s1, v2 = e2-s2, d = s2-s1; - auto a = v1.cross(v2), a1 = v1.cross(d), a2 = v2.cross(d); - if (a == 0) { //if parallel - auto b1=s1.dot(v1), c1=e1.dot(v1), - b2=s2.dot(v1), c2=e2.dot(v1); - if (a1 || a2 || max(b1,min(b2,c2))>min(c1,max(b2,c2))) - return 0; - r1 = min(b2,c2)c1 ? e1 : (b2>c2 ? s2 : e2); - return 2-(r1==r2); - } - if (a < 0) { a = -a; a1 = -a1; a2 = -a2; } - if (0min(c1,max(b2,c2))) + return 0; + r1 = min(b2,c2)c1 ? e1 : (b2>c2 ? s2 : e2); + return 2-(r1==r2); + } + if (a < 0) { a = -a; a1 = -a1; a2 = -a2; } + if (0 P; bool eq(P a, P b) { - return (a-b).dist()<1e-8; + return (a-b).dist()<1e-8; } int main() { - rep(t,0,1000000) { - const int GRID=6; - P a(rand()%GRID, rand()%GRID), b(rand()%GRID, rand()%GRID), c(rand()%GRID, rand()%GRID), d(rand()%GRID, rand()%GRID); - P tmp1, tmp2; - auto res = oldImpl::segmentIntersection(a,b,c,d, tmp1, tmp2); - auto res2 = segInter(a,b,c,d); - if (res != sz(res2)) { - cout< a(res2.begin(), res2.end()); - vector

b({tmp1, tmp2}); - sort(all(b)); - assert(eq(a[0], b[0]) && eq(a[1],b[1])); - } - } - cout<<"Tests passed!"< a(res2.begin(), res2.end()); + vector

b({tmp1, tmp2}); + sort(all(b)); + assert(eq(a[0], b[0]) && eq(a[1],b[1])); + } + } + cout<<"Tests passed!"< P; bool eq(P a, P b) { - return (a-b).dist() poly; - rep(j,0, numPts) - poly.push_back(P(rand()%range, rand()%range)); - poly = genPolygon(poly); - rep(i,0,PTPERPOLY){ - P p(rand()%range, rand()%range); - assert(inPolygon(poly, p, true) == old::insidePolygon(all(poly), p, true)); - assert(inPolygon(poly, p, false) == old::insidePolygon(all(poly), p, false)); - } - } + rep(i,0,NUMPOLY) { + vector

poly; + rep(j,0, numPts) + poly.push_back(P(rand()%range, rand()%range)); + poly = genPolygon(poly); + rep(i,0,PTPERPOLY){ + P p(rand()%range, rand()%range); + assert(inPolygon(poly, p, true) == old::insidePolygon(all(poly), p, true)); + assert(inPolygon(poly, p, false) == old::insidePolygon(all(poly), p, false)); + } + } } int main() { - test(20,5); - test(1001,100); - test(1000,1000); - cout<<"Tests passed!"<> edges; - vector mat(n, vi(n)); - rep(it,0,m) { - int i = rand() % n; - int j = rand() % n; - if (i == j) continue; - int w = rand() % mxFlow; - edges.push_back({i, j, w}); - mat[i][j] += w; - mat[j][i] += w; - } - auto calc = [&](int s, int t) { - Dinic flow(n); - for (auto e : edges) { - flow.addEdge((int)e[0], (int)e[1], e[2], e[2]); - } - return flow.calc(s, t); - }; - vector gomoryHuTree = gomoryHu(n, edges); - vector>> adj(n); - for (auto e : gomoryHuTree) { - adj[e[0]].push_back({(int)e[1], (int)e[2]}); - adj[e[1]].push_back({(int)e[0], (int)e[2]}); - } - auto dfs = make_y_combinator([&](auto dfs, int start, int cur, int p, int mn) -> void { - if (start != cur) { - assert(mn == calc(start, cur)); - } - for (auto i : adj[cur]) { - if (i[0] != p) - dfs(start, i[0], cur, min(mn, i[1])); - } - }); - dfs(0, 0, -1, INT_MAX); + for (int it = 0; it < iters; it++) { + int n = rand()%N+1; + int m = rand()%(N*N); + vector> edges; + vector mat(n, vi(n)); + rep(it,0,m) { + int i = rand() % n; + int j = rand() % n; + if (i == j) continue; + int w = rand() % mxFlow; + edges.push_back({i, j, w}); + mat[i][j] += w; + mat[j][i] += w; + } + auto calc = [&](int s, int t) { + Dinic flow(n); + for (auto e : edges) { + flow.addEdge((int)e[0], (int)e[1], e[2], e[2]); + } + return flow.calc(s, t); + }; + vector gomoryHuTree = gomoryHu(n, edges); + vector>> adj(n); + for (auto e : gomoryHuTree) { + adj[e[0]].push_back({(int)e[1], (int)e[2]}); + adj[e[1]].push_back({(int)e[0], (int)e[2]}); + } + auto dfs = make_y_combinator([&](auto dfs, int start, int cur, int p, int mn) -> void { + if (start != cur) { + assert(mn == calc(start, cur)); + } + for (auto i : adj[cur]) { + if (i[0] != p) + dfs(start, i[0], cur, min(mn, i[1])); + } + }); + dfs(0, 0, -1, INT_MAX); - // Check that the lightest edge agrees with GlobalMinCut. - if (n >= 2) { - ll minCut = LLONG_MAX; - for (auto e : gomoryHuTree) { - minCut = min(minCut, e[2]); - } - auto mat2 = mat; - auto pa = globalMinCut(mat2); - assert(pa.first == minCut); - vi inCut(n); - assert(sz(pa.second) != 0); - assert(sz(pa.second) != n); - for (int x : pa.second) { - assert(0 <= x && x < n); - assert(!inCut[x]); - inCut[x] = 1; - } - int cutw = 0; - rep(i,0,n) rep(j,0,n) if (inCut[i] && !inCut[j]) { - cutw += mat[i][j]; - } - assert(pa.first == cutw); - } - } + // Check that the lightest edge agrees with GlobalMinCut. + if (n >= 2) { + ll minCut = LLONG_MAX; + for (auto e : gomoryHuTree) { + minCut = min(minCut, e[2]); + } + auto mat2 = mat; + auto pa = globalMinCut(mat2); + assert(pa.first == minCut); + vi inCut(n); + assert(sz(pa.second) != 0); + assert(sz(pa.second) != n); + for (int x : pa.second) { + assert(0 <= x && x < n); + assert(!inCut[x]); + inCut[x] = 1; + } + int cutw = 0; + rep(i,0,n) rep(j,0,n) if (inCut[i] && !inCut[j]) { + cutw += mat[i][j]; + } + assert(pa.first == cutw); + } + } } signed main() { - test(25, 5, 200); - test(100, 1000, 5); - test(100, 1, 20); - test(5, 5, 20000); - cout<<"Tests passed!"<> tree; - vector vals; - vector pars; - int unit = -1e9; - int f(int a, int b) { return max(a, b); } - void root(int cur, int p = -1) { - pars[cur] = p; - for (auto i: tree[cur]) { - if (i != p) root(i, cur); - } - } - bruteforce(vector> _tree): tree(_tree), vals(sz(tree)), pars(sz(tree)) { - root(0); - } - bool dfsModify(int cur, int target, int val, int p=-1) { - if (cur == target) { - vals[cur] += val; - return true; - } - bool alongPath = false; - for (auto i: tree[cur]) { - if (i == p) continue; - alongPath |= dfsModify(i, target, val, cur); - } - if (alongPath) vals[cur] += val; - return alongPath; - } - void modifyPath(int a, int b, int val) { - dfsModify(a, b, val); - } + vector> tree; + vector vals; + vector pars; + int unit = -1e9; + int f(int a, int b) { return max(a, b); } + void root(int cur, int p = -1) { + pars[cur] = p; + for (auto i: tree[cur]) { + if (i != p) root(i, cur); + } + } + bruteforce(vector> _tree): tree(_tree), vals(sz(tree)), pars(sz(tree)) { + root(0); + } + bool dfsModify(int cur, int target, int val, int p=-1) { + if (cur == target) { + vals[cur] += val; + return true; + } + bool alongPath = false; + for (auto i: tree[cur]) { + if (i == p) continue; + alongPath |= dfsModify(i, target, val, cur); + } + if (alongPath) vals[cur] += val; + return alongPath; + } + void modifyPath(int a, int b, int val) { + dfsModify(a, b, val); + } - int dfsQuery(int cur, int target, int p = -1) { - if (cur == target) { - return vals[cur]; - } - int res = unit; - for (auto i: tree[cur]) { - if (i == p) continue; - res = f(res, dfsQuery(i, target, cur)); - } - if (res != unit) { - return f(res, vals[cur]); - } - return res; - } - int queryPath(int a, int b) { - return dfsQuery(a, b); - } - int dfsSubtree(int cur, int p) { - int res = vals[cur]; - for (auto i: tree[cur]) { - if (i != p) - res = f(res, dfsSubtree(i, cur)); - } - return res; - } - int querySubtree(int a) { - return dfsSubtree(a, pars[a]); - } + int dfsQuery(int cur, int target, int p = -1) { + if (cur == target) { + return vals[cur]; + } + int res = unit; + for (auto i: tree[cur]) { + if (i == p) continue; + res = f(res, dfsQuery(i, target, cur)); + } + if (res != unit) { + return f(res, vals[cur]); + } + return res; + } + int queryPath(int a, int b) { + return dfsQuery(a, b); + } + int dfsSubtree(int cur, int p) { + int res = vals[cur]; + for (auto i: tree[cur]) { + if (i != p) + res = f(res, dfsSubtree(i, cur)); + } + return res; + } + int querySubtree(int a) { + return dfsSubtree(a, pars[a]); + } }; void testAgainstOld(int n, int iters, int queries) { - for (int trees = 0; trees < iters; trees++) { - auto graph = genRandomTree(n); - vector> tree1(n); - vector>> tree2(n); - for (auto i : graph) { - tree1[i.first].push_back(i.second); - tree1[i.second].push_back(i.first); - } - for (int i = 0; i < sz(tree1); i++) { - for (auto j : tree1[i]) { - tree2[i].push_back({j, 0}); - } - } - HLD hld(tree1); - old::HLD hld2(tree2); - hld.tree->set(0, n, 0); - for (int itr = 0; itr < queries; itr++) { - if (rand() % 2) { - int node = rand() % n; - int val = rand() % 10; - hld2.update(node, val); - hld.modifyPath(node, node, val - hld.queryPath(node, node)); - } else { - int a = rand() % n; - int b = rand() % n; - assert(hld.queryPath(a, b) == hld2.query2(a, b).first); - } - } - } + for (int trees = 0; trees < iters; trees++) { + auto graph = genRandomTree(n); + vector> tree1(n); + vector>> tree2(n); + for (auto i : graph) { + tree1[i.first].push_back(i.second); + tree1[i.second].push_back(i.first); + } + for (int i = 0; i < sz(tree1); i++) { + for (auto j : tree1[i]) { + tree2[i].push_back({j, 0}); + } + } + HLD hld(tree1); + old::HLD hld2(tree2); + hld.tree->set(0, n, 0); + for (int itr = 0; itr < queries; itr++) { + if (rand() % 2) { + int node = rand() % n; + int val = rand() % 10; + hld2.update(node, val); + hld.modifyPath(node, node, val - hld.queryPath(node, node)); + } else { + int a = rand() % n; + int b = rand() % n; + assert(hld.queryPath(a, b) == hld2.query2(a, b).first); + } + } + } } void testAgainstBrute(int n, int iters, int queries) { - for (int trees = 0; trees < iters; trees++) { - auto graph = genRandomTree(n); - vector> tree1(n); - for (auto i : graph) { - tree1[i.first].push_back(i.second); - tree1[i.second].push_back(i.first); - } - HLD hld(tree1); - bruteforce hld2(tree1); - hld.tree->set(0, n, 0); - for (int itr = 0; itr < queries; itr++) { - int rng = rand() % 3; - if (rng == 0) { - int a = rand() % n; - int b = rand() % n; - int val = rand() % 10; - hld.modifyPath(a, b, val); - hld2.modifyPath(a, b, val); - } else if (rng == 1){ - int a = rand() % n; - int b = rand() % n; - hld.queryPath(a, b); - hld2.queryPath(a, b); - assert(hld.queryPath(a, b) == hld2.queryPath(a, b)); - } else if (rng == 2) { - int a = rand() % n; - assert(hld.querySubtree(a) == hld2.querySubtree(a)); - } - } - } + for (int trees = 0; trees < iters; trees++) { + auto graph = genRandomTree(n); + vector> tree1(n); + for (auto i : graph) { + tree1[i.first].push_back(i.second); + tree1[i.second].push_back(i.first); + } + HLD hld(tree1); + bruteforce hld2(tree1); + hld.tree->set(0, n, 0); + for (int itr = 0; itr < queries; itr++) { + int rng = rand() % 3; + if (rng == 0) { + int a = rand() % n; + int b = rand() % n; + int val = rand() % 10; + hld.modifyPath(a, b, val); + hld2.modifyPath(a, b, val); + } else if (rng == 1){ + int a = rand() % n; + int b = rand() % n; + hld.queryPath(a, b); + hld2.queryPath(a, b); + assert(hld.queryPath(a, b) == hld2.queryPath(a, b)); + } else if (rng == 2) { + int a = rand() % n; + assert(hld.querySubtree(a) == hld2.querySubtree(a)); + } + } + } } int main() { - srand(2); - testAgainstBrute(5, 1000, 10000); - testAgainstBrute(1000, 100, 100); - testAgainstOld(5, 1000, 100); - testAgainstOld(10000, 100, 1000); - cout<<"Tests passed!"< vpi; typedef vector graph; struct LCA { - vi time; - vector dist; - RMQ rmq; + vi time; + vector dist; + RMQ rmq; - LCA(graph& C) : time(sz(C), -99), dist(sz(C)), rmq(dfs(C)) {} + LCA(graph& C) : time(sz(C), -99), dist(sz(C)), rmq(dfs(C)) {} - vpi dfs(graph& C) { - vector> q(1); - vpi ret; - int T = 0, v, p, d; ll di; - while (!q.empty()) { - tie(v, p, d, di) = q.back(); - q.pop_back(); - if (d) ret.emplace_back(d, p); - time[v] = T++; - dist[v] = di; - for(auto &e: C[v]) if (e.first != p) - q.emplace_back(e.first, v, d+1, di + e.second); - } - return ret; - } + vpi dfs(graph& C) { + vector> q(1); + vpi ret; + int T = 0, v, p, d; ll di; + while (!q.empty()) { + tie(v, p, d, di) = q.back(); + q.pop_back(); + if (d) ret.emplace_back(d, p); + time[v] = T++; + dist[v] = di; + for(auto &e: C[v]) if (e.first != p) + q.emplace_back(e.first, v, d+1, di + e.second); + } + return ret; + } - int query(int a, int b) { - if (a == b) return a; - a = time[a], b = time[b]; - return rmq.query(min(a, b), max(a, b)).second; - } - ll distance(int a, int b) { - int lca = query(a, b); - return dist[a] + dist[b] - 2 * dist[lca]; - } + int query(int a, int b) { + if (a == b) return a; + a = time[a], b = time[b]; + return rmq.query(min(a, b), max(a, b)).second; + } + ll distance(int a, int b) { + int lca = query(a, b); + return dist[a] + dist[b] - 2 * dist[lca]; + } }; } void getPars(vector &tree, int cur, int p, int d, vector &par, vector &depth) { - par[cur] = p; - depth[cur] = d; - for(auto i: tree[cur]) if (i != p) { - getPars(tree, i, cur, d+1, par, depth); - } + par[cur] = p; + depth[cur] = d; + for(auto i: tree[cur]) if (i != p) { + getPars(tree, i, cur, d+1, par, depth); + } } void test_n(int n, int num) { - for (int out=0; out tree(n); - vector>> oldTree(n); - for (auto i: graph) { - tree[i.first].push_back(i.second); - tree[i.second].push_back(i.first); - oldTree[i.first].push_back({i.second, 1}); - oldTree[i.second].push_back({i.first, 1}); - } - vector par(n), depth(n); - getPars(tree, 0, 0, 0, par, depth); - vector tbl = treeJump(par); - LCA new_lca(tree); - old::LCA old_lca(oldTree); - for (int i=0; i<100; i++) { - int a = rand()%n, b = rand()%n; - int binLca = lca(tbl, depth, a, b); - int newLca = new_lca.lca(a,b); - int oldLca = old_lca.query(a,b); - assert(oldLca == newLca); - assert(binLca == newLca); - } - } + for (int out=0; out tree(n); + vector>> oldTree(n); + for (auto i: graph) { + tree[i.first].push_back(i.second); + tree[i.second].push_back(i.first); + oldTree[i.first].push_back({i.second, 1}); + oldTree[i.second].push_back({i.first, 1}); + } + vector par(n), depth(n); + getPars(tree, 0, 0, 0, par, depth); + vector tbl = treeJump(par); + LCA new_lca(tree); + old::LCA old_lca(oldTree); + for (int i=0; i<100; i++) { + int a = rand()%n, b = rand()%n; + int binLca = lca(tbl, depth, a, b); + int newLca = new_lca.lca(a,b); + int oldLca = old_lca.query(a,b); + assert(oldLca == newLca); + assert(binLca == newLca); + } + } } signed main() { - test_n(10, 1000); - test_n(100, 100); - test_n(1000, 10); - cout<<"Tests passed!"< ed2(n); - int p =rand()%100; - rep(i, 0, n) rep(j, 0, i) { - ed[i][j] = (rand() % 100) < p; - ed[j][i] = ed[i][j]; - ed2[i][j] = ed[i][j]; - ed2[j][i] = ed[j][i]; - } - Maxclique clique2(ed); - int mx = 0; - maximal::cliques(ed2, [&](auto x){mx = max(mx, int(x.count()));}); - assert(mx == sz(clique2.maxClique())); - } - cout<<"Tests passed!"< ed2(n); + int p =rand()%100; + rep(i, 0, n) rep(j, 0, i) { + ed[i][j] = (rand() % 100) < p; + ed[j][i] = ed[i][j]; + ed2[i][j] = ed[i][j]; + ed2[j][i] = ed[j][i]; + } + Maxclique clique2(ed); + int mx = 0; + maximal::cliques(ed2, [&](auto x){mx = max(mx, int(x.count()));}); + assert(mx == sz(clique2.maxClique())); + } + cout<<"Tests passed!"< m) - swap(n, m); + for (int it = 0; it < iters; it++) { + int n = randRange(0, N), m = randRange(0, N); + if (n > m) + swap(n, m); - MCMF mcmf(n + m + 2); - int s = 0; - int t = 1; - for (int i = 0; i < n; i++) - mcmf.addEdge(s, i + 2, 1, 0); - for (int i = 0; i < m; i++) - mcmf.addEdge(2 + n + i, t, 1, 0); + MCMF mcmf(n + m + 2); + int s = 0; + int t = 1; + for (int i = 0; i < n; i++) + mcmf.addEdge(s, i + 2, 1, 0); + for (int i = 0; i < m; i++) + mcmf.addEdge(2 + n + i, t, 1, 0); - vector cost(n, vi(m)); - for (int i = 0; i < n; i++) { - for (int j = 0; j < m; j++) { - cost[i][j] = randRange(-mxCost, mxCost); - mcmf.addEdge(i + 2, 2 + n + j, 1, cost[i][j]); - } - } - mcmf.setpi(s); - auto maxflow = mcmf.maxflow(s, t); - auto matching = hungarian(cost); - assert(maxflow.first == n); - assert(maxflow.second == matching.first); - int matchSum = 0; - set used; - for (int i = 0; i < n; i++) { - matchSum += cost[i][matching.second[i]]; - assert(used.count(matching.second[i]) == 0); - used.insert(matching.second[i]); - } - assert(matchSum == matching.first); - return; - } + vector cost(n, vi(m)); + for (int i = 0; i < n; i++) { + for (int j = 0; j < m; j++) { + cost[i][j] = randRange(-mxCost, mxCost); + mcmf.addEdge(i + 2, 2 + n + j, 1, cost[i][j]); + } + } + mcmf.setpi(s); + auto maxflow = mcmf.maxflow(s, t); + auto matching = hungarian(cost); + assert(maxflow.first == n); + assert(maxflow.second == matching.first); + int matchSum = 0; + set used; + for (int i = 0; i < n; i++) { + matchSum += cost[i][matching.second[i]]; + assert(used.count(matching.second[i]) == 0); + used.insert(matching.second[i]); + } + assert(matchSum == matching.first); + return; + } } signed main() { - test(25, 5, 1000); - test(100, 1000, 100); - test(100, 1, 50); - test(5, 5, 10000); - cout << "Tests passed!" << endl; + test(25, 5, 1000); + test(100, 1000, 100); + test(100, 1, 50); + test(5, 5, 10000); + cout << "Tests passed!" << endl; } diff --git a/stress-tests/number-theory/Factor.cpp b/stress-tests/number-theory/Factor.cpp index 0ac3cdd22..d5b982d38 100644 --- a/stress-tests/number-theory/Factor.cpp +++ b/stress-tests/number-theory/Factor.cpp @@ -4,33 +4,33 @@ mt19937_64 uni(time(0)); void assertValid(ull N, vector prFac){ - ull cur=1; - for (auto i: prFac){ - if (!isPrime(i)){ - cout<{2}); - assert((factor(2299) == vector{11, 19, 11})); - rep(n,2,1e5) { - auto res = factor(n); - assertValid(n, res); - res = factor(n*ll(n)); - assertValid(n*ll(n), res); - } - rep(i,2,1e5) { - ull n = 1 + (uni()%(3ul<<61)); - auto res = factor(n); - assertValid(n, res); - } - cout<<"Tests passed!"<{2}); + assert((factor(2299) == vector{11, 19, 11})); + rep(n,2,1e5) { + auto res = factor(n); + assertValid(n, res); + res = factor(n*ll(n)); + assertValid(n*ll(n), res); + } + rep(i,2,1e5) { + ull n = 1 + (uni()%(3ul<<61)); + auto res = factor(n); + assertValid(n, res); + } + cout<<"Tests passed!"<> 1); - return e & 1 ? x * a % mod : x; + if (e == 0) + return 1; + ll x = modpow(a * a % mod, e >> 1); + return e & 1 ? x * a % mod : x; } vl simpleConv(vl a, vl b) { @@ -24,11 +24,11 @@ vl simpleConv(vl a, vl b) { } int ra() { - static unsigned X; - X *= 123671231; - X += 1238713; - X ^= 1237618; - return (X >> 1); + static unsigned X; + X *= 123671231; + X += 1238713; + X ^= 1237618; + return (X >> 1); } int main() { @@ -43,13 +43,13 @@ int main() { for(auto &x: simpleConv(a, b)) res += (ll)x * ind++ % mod; for(auto &x: conv(a, b)) res2 += (ll)x * ind2++ % mod; a.resize(16); - vl a2 = a; - ntt(a2); - rep(k, 0, sz(a2)) { - ll sum = 0; - rep(x, 0, sz(a2)) { sum = (sum + a[x] * modpow(root, k * x * (mod - 1) / sz(a))) % mod; } - assert(sum == a2[k]); - } + vl a2 = a; + ntt(a2); + rep(k, 0, sz(a2)) { + ll sum = 0; + rep(x, 0, sz(a2)) { sum = (sum + a[x] * modpow(root, k * x * (mod - 1) / sz(a))) % mod; } + assert(sum == a2[k]); + } } assert(res==res2); cout<<"Tests passed!"< void gen(string &s, int at, int alpha, F f) { - if (at == sz(s)) - f(); - else { - rep(i, 0, alpha) { - s[at] = (char)('a' + i); - gen(s, at + 1, alpha, f); - } - } + if (at == sz(s)) + f(); + else { + rep(i, 0, alpha) { + s[at] = (char)('a' + i); + gen(s, at + 1, alpha, f); + } + } } void test(const string &s) { - int n = sz(s); - vi found = Z(s); - vi expected(n, 0); - rep(i, 1, n) { // exclude index 0 (!) - int j = 0; - while (i + j < n && s[i + j] == s[j]) - j++; - expected[i] = j; - } - assert(found == expected); + int n = sz(s); + vi found = Z(s); + vi expected(n, 0); + rep(i, 1, n) { // exclude index 0 (!) + int j = 0; + while (i + j < n && s[i + j] == s[j]) + j++; + expected[i] = j; + } + assert(found == expected); } signed main() { - ios::sync_with_stdio(0); - cin.tie(0); - rep(n, 0, 13) { - string s(n, 'x'); - gen(s, 0, 3, [&]() { test(s); }); - } - rep(n, 0, 11) { - string s(n, 'x'); - gen(s, 0, 4, [&]() { test(s); }); - } - cout<<"Tests passed!"<(end - begin).count(); - cerr << duration << "ms elapsed [" << label << "]" << endl; - } + decltype(chrono::high_resolution_clock::now()) begin; + const string label; + timeit(string label = "???") : label(label) { begin = chrono::high_resolution_clock::now(); } + ~timeit() { + auto end = chrono::high_resolution_clock::now(); + auto duration = chrono::duration_cast(end - begin).count(); + cerr << duration << "ms elapsed [" << label << "]" << endl; + } }; diff --git a/stress-tests/utilities/genPolygon.h b/stress-tests/utilities/genPolygon.h index ee81c93f3..641f804e9 100644 --- a/stress-tests/utilities/genPolygon.h +++ b/stress-tests/utilities/genPolygon.h @@ -7,66 +7,66 @@ #include "random.h" template pair> conquer(vector

pts, int depth) { - if (depth>100) { - return {false, {}}; - } - if (sz(pts) <= 2) return {true,pts}; - if (sz(pts) == 3) { - swap(pts[1], pts[2]); - return {true,pts}; - } + if (depth>100) { + return {false, {}}; + } + if (sz(pts) <= 2) return {true,pts}; + if (sz(pts) == 3) { + swap(pts[1], pts[2]); + return {true,pts}; + } - int divideId = randRange(2, sz(pts)); - P p1 = pts[divideId]; - double divideK = randDouble(0.01, 0.99); - P p2(divideK*(pts[1].x-pts[0].x) + pts[0].x, divideK*(pts[1].y - pts[0].y) + pts[0].y); - vector line = {p2.y - p1.y, p1.x - p2.x, -p1.x*p2.y + p1.y*p2.x}; - int idx0 = ((line[0]*pts[0].x + line[1]*pts[0].y + line[2]) >=0); - int idx1 = ((line[0]*pts[1].x + line[1]*pts[1].y + line[2]) >=0); - if (idx0==idx1) - return conquer(pts, depth+1); - array, 2> S; - S[idx0].push_back(pts[0]); - S[idx0].push_back(p1); - S[!idx0].push_back(p1); - S[!idx0].push_back(pts[1]); - rep(i,2,sz(pts)) { - if (i == divideId) continue; - int idx = ((line[0]*pts[i].x + line[1]*pts[i].y + line[2]) >=0); - S[idx].push_back(pts[i]); - } - auto pa = conquer(S[idx0], depth+1); - auto pb = conquer(S[!idx0], depth+1); - if (!pa.first || !pb.first) return {false, {}}; - pb.second.erase(pb.second.begin()); - pa.second.insert(pa.second.end(), all(pb.second)); - return pa; + int divideId = randRange(2, sz(pts)); + P p1 = pts[divideId]; + double divideK = randDouble(0.01, 0.99); + P p2(divideK*(pts[1].x-pts[0].x) + pts[0].x, divideK*(pts[1].y - pts[0].y) + pts[0].y); + vector line = {p2.y - p1.y, p1.x - p2.x, -p1.x*p2.y + p1.y*p2.x}; + int idx0 = ((line[0]*pts[0].x + line[1]*pts[0].y + line[2]) >=0); + int idx1 = ((line[0]*pts[1].x + line[1]*pts[1].y + line[2]) >=0); + if (idx0==idx1) + return conquer(pts, depth+1); + array, 2> S; + S[idx0].push_back(pts[0]); + S[idx0].push_back(p1); + S[!idx0].push_back(p1); + S[!idx0].push_back(pts[1]); + rep(i,2,sz(pts)) { + if (i == divideId) continue; + int idx = ((line[0]*pts[i].x + line[1]*pts[i].y + line[2]) >=0); + S[idx].push_back(pts[i]); + } + auto pa = conquer(S[idx0], depth+1); + auto pb = conquer(S[!idx0], depth+1); + if (!pa.first || !pb.first) return {false, {}}; + pb.second.erase(pb.second.begin()); + pa.second.insert(pa.second.end(), all(pb.second)); + return pa; } template vector

genPolygon(vector

pts, int depth=0) { - if (depth>100) return {P(0,0), P(1,0), P(0,1)}; - sort(all(pts)); - pts.resize(unique(all(pts)) - pts.begin()); - shuffle_vec(pts); - if (sz(pts) <=3) return pts; - vector line ={(double)(pts[1].y-pts[0].y), (double)(pts[0].x - pts[1].x), (double)(-pts[0].x*pts[1].y + pts[0].y*pts[1].x)}; - array, 2> S; - S[0].push_back(pts[0]); - S[0].push_back(pts[1]); - S[1].push_back(pts[1]); - S[1].push_back(pts[0]); - rep(i,2,sz(pts)) { - int idx = (line[0]*pts[i].x + line[1]*pts[i].y + line[2]) >=0; - S[idx].push_back(pts[i]); - } - auto ta = conquer(S[0],0); - auto tb = conquer(S[1],0); - auto pa=ta.second, pb=tb.second; - if (!ta.first || !tb.first) { - return genPolygon(pts, depth+1); - } - pa.erase(pa.begin()); - pb.erase(pb.begin()); - pa.insert(pa.end(), all(pb)); - if (polygonArea2(pa) < 0) reverse(all(pa)); - return pa; + if (depth>100) return {P(0,0), P(1,0), P(0,1)}; + sort(all(pts)); + pts.resize(unique(all(pts)) - pts.begin()); + shuffle_vec(pts); + if (sz(pts) <=3) return pts; + vector line ={(double)(pts[1].y-pts[0].y), (double)(pts[0].x - pts[1].x), (double)(-pts[0].x*pts[1].y + pts[0].y*pts[1].x)}; + array, 2> S; + S[0].push_back(pts[0]); + S[0].push_back(pts[1]); + S[1].push_back(pts[1]); + S[1].push_back(pts[0]); + rep(i,2,sz(pts)) { + int idx = (line[0]*pts[i].x + line[1]*pts[i].y + line[2]) >=0; + S[idx].push_back(pts[i]); + } + auto ta = conquer(S[0],0); + auto tb = conquer(S[1],0); + auto pa=ta.second, pb=tb.second; + if (!ta.first || !tb.first) { + return genPolygon(pts, depth+1); + } + pa.erase(pa.begin()); + pb.erase(pb.begin()); + pa.insert(pa.end(), all(pb)); + if (polygonArea2(pa) < 0) reverse(all(pa)); + return pa; } diff --git a/stress-tests/utilities/genTree.h b/stress-tests/utilities/genTree.h index 1cb6939f5..97d5a312a 100644 --- a/stress-tests/utilities/genTree.h +++ b/stress-tests/utilities/genTree.h @@ -6,53 +6,53 @@ */ vector> pruferCodeToTree(vector &pruferCode) { - // Stores number count of nodes in the prufer code - unordered_map nodeCount; + // Stores number count of nodes in the prufer code + unordered_map nodeCount; - // Set of integers absent in prufer code. They are the leaves - set leaves; + // Set of integers absent in prufer code. They are the leaves + set leaves; - int len = (int) pruferCode.size(); - int node = len + 2; + int len = (int) pruferCode.size(); + int node = len + 2; - // Count frequency of nodes - for ( int i = 0; i < len; i++ ) { - int t = pruferCode[i]; - nodeCount[t]++; - } + // Count frequency of nodes + for ( int i = 0; i < len; i++ ) { + int t = pruferCode[i]; + nodeCount[t]++; + } - // Find the absent nodes - for ( int i = 1; i <= node; i++ ) { - if ( nodeCount.find ( i ) == nodeCount.end() ) leaves.insert ( i ); - } + // Find the absent nodes + for ( int i = 1; i <= node; i++ ) { + if ( nodeCount.find ( i ) == nodeCount.end() ) leaves.insert ( i ); + } - vector> edges; - /*Connect Edges*/ - for ( int i = 0; i < len; i++ ){ - int a = pruferCode[i]; // First node + vector> edges; + /*Connect Edges*/ + for ( int i = 0; i < len; i++ ){ + int a = pruferCode[i]; // First node - //Find the smallest number which is not present in prufer code now - int b = *leaves.begin(); // the leaf + //Find the smallest number which is not present in prufer code now + int b = *leaves.begin(); // the leaf - edges.push_back({a,b}); // Edge of the tree + edges.push_back({a,b}); // Edge of the tree - leaves.erase ( b ); // Remove from absent list - nodeCount[a]--; // Remove from prufer code - if ( nodeCount[a] == 0 ) leaves.insert ( a ); // If a becomes absent - } + leaves.erase ( b ); // Remove from absent list + nodeCount[a]--; // Remove from prufer code + if ( nodeCount[a] == 0 ) leaves.insert ( a ); // If a becomes absent + } - // The final edge - edges.push_back({*leaves.begin(), *leaves.rbegin()}); - return edges; + // The final edge + edges.push_back({*leaves.begin(), *leaves.rbegin()}); + return edges; } vector> genRandomTree(int n) { - vector pruferCode; - for (int i=0; i pruferCode; + for (int i=0; i Point randIntPt(int lim) { - return Point(rand()%(lim*2) - lim, rand()%(lim*2)-lim); + return Point(rand()%(lim*2) - lim, rand()%(lim*2)-lim); } diff --git a/stress-tests/utilities/random.h b/stress-tests/utilities/random.h index 7b4f0a7c8..196914057 100644 --- a/stress-tests/utilities/random.h +++ b/stress-tests/utilities/random.h @@ -4,7 +4,7 @@ // returns random int in [0, hi), like Python's random.randrange int randRange(int hi) { - return rand() % hi; + return rand() % hi; } bool randBool() { @@ -13,109 +13,109 @@ bool randBool() { // returns random int in [lo, hi), like Python's random.randrange int randRange(int lo, int hi) { - return lo + randRange(hi - lo); + return lo + randRange(hi - lo); } // returns random int in [0, 2^64) uint64_t randU64() { - uint64_t a = rand() & 0xffffff; - uint64_t b = rand() & 0xffffff; - uint64_t c = rand() & 0xffff; - return a << 40 | b << 16 | c; + uint64_t a = rand() & 0xffffff; + uint64_t b = rand() & 0xffffff; + uint64_t c = rand() & 0xffff; + return a << 40 | b << 16 | c; } uint64_t randRange(uint64_t hi) { - return randU64() % hi; + return randU64() % hi; } uint64_t randRange(uint64_t lo, uint64_t hi) { - return lo + randRange(hi - lo); + return lo + randRange(hi - lo); } int64_t randRange(int64_t hi) { - assert(hi > 0); - return (int64_t) randRange((uint64_t) hi); + assert(hi > 0); + return (int64_t) randRange((uint64_t) hi); } int64_t randRange(int64_t lo, int64_t hi) { - return lo + randRange(hi - lo); + return lo + randRange(hi - lo); } // returns random int in [0, hi], like Python's random.randint int randIncl(int hi) { - return randRange(hi + 1); + return randRange(hi + 1); } // returns random int in [lo, hi], like Python's random.randint int randIncl(int lo, int hi) { - return randRange(lo, hi + 1); + return randRange(lo, hi + 1); } int64_t randIncl(int64_t hi) { - return randRange(hi + 1); + return randRange(hi + 1); } int64_t randIncl(int64_t lo, int64_t hi) { - return randRange(lo, hi + 1); + return randRange(lo, hi + 1); } // returns uniformly random double in [lo, hi) double randDouble(double lo, double hi) { - static mt19937_64 rng(time(0)); - std::uniform_real_distribution<> dis(lo, hi); - return dis(rng); + static mt19937_64 rng(time(0)); + std::uniform_real_distribution<> dis(lo, hi); + return dis(rng); } // int -> double based on IEEE 754 bitpattern double bitPatternToDouble(uint64_t x) { - union { - double d; - uint64_t i; - } u; - u.i = x; - return u.d; + union { + double d; + uint64_t i; + } u; + u.i = x; + return u.d; } // double -> int based on IEEE 754 bitpattern uint64_t doubleToBitPattern(double x) { - union { - double d; - uint64_t i; - } u; - u.d = x; - return u.i; + union { + double d; + uint64_t i; + } u; + u.d = x; + return u.i; } // random double double randDoubleUniformBitPattern() { - return bitPatternToDouble(randU64()); + return bitPatternToDouble(randU64()); } // random double in [lo, hi), with any bit pattern being equally likely double randDoubleUniformBitPattern(double lo, double hi) { - return bitPatternToDouble(randRange(doubleToBitPattern(lo), doubleToBitPattern(hi))); + return bitPatternToDouble(randRange(doubleToBitPattern(lo), doubleToBitPattern(hi))); } // add ~y ulps (units of last precision) to x, similar to calling next_after y times double addUlps(double x, int64_t y) { - if (x == 0 && y < 0) { - return -addUlps(-x, -y); - } - return bitPatternToDouble(doubleToBitPattern(x) + y); + if (x == 0 && y < 0) { + return -addUlps(-x, -y); + } + return bitPatternToDouble(doubleToBitPattern(x) + y); } // random int in [-lim, lim], perturbed by a few ulps double randNearIntUlps(int lim, int64_t ulps = 5) { - return addUlps(randIncl(-lim, lim), randIncl(-ulps, ulps)); + return addUlps(randIncl(-lim, lim), randIncl(-ulps, ulps)); } // random int in [-lim, lim], perturbed by a random double in [-eps, eps] double randNearIntEps(int lim, double eps) { - return randIncl(-lim, lim) + randDouble(-eps, eps); + return randIncl(-lim, lim) + randDouble(-eps, eps); } // like random_shuffle but uses rand() as RNG source template void shuffle_vec(T& vec) { - random_shuffle(begin(vec), end(vec), [](int lim) { return rand() % lim; }); + random_shuffle(begin(vec), end(vec), [](int lim) { return rand() % lim; }); } diff --git a/stress-tests/utilities/utils.h b/stress-tests/utilities/utils.h index 7da69501a..1c728e01a 100644 --- a/stress-tests/utilities/utils.h +++ b/stress-tests/utilities/utils.h @@ -2,13 +2,13 @@ // Taken from https://stackoverflow.com/a/40873657/2014771 template struct y_combinator { - F f; // the lambda will be stored here + F f; // the lambda will be stored here - // a forwarding operator(): - template decltype(auto) operator()(Args &&... args) const { - // we pass ourselves to f, then the arguments. - return f(std::ref(*this), std::forward(args)...); - } + // a forwarding operator(): + template decltype(auto) operator()(Args &&... args) const { + // we pass ourselves to f, then the arguments. + return f(std::ref(*this), std::forward(args)...); + } }; // helper function that deduces the type of the lambda: diff --git a/stress-tests/various/FastKnapsack.cpp b/stress-tests/various/FastKnapsack.cpp index 508d4b4c6..289a25e85 100644 --- a/stress-tests/various/FastKnapsack.cpp +++ b/stress-tests/various/FastKnapsack.cpp @@ -3,30 +3,30 @@ #include "../../content/various/FastKnapsack.h" int naive(vi w, int t) { - vector can_reach(t+1); - can_reach[0] = true; - for (int x : w) { - for (int i = t-x; i >= 0; --i) { - if (can_reach[i]) can_reach[i+x] = true; - } - } - for (int i = t;; i--) - if (can_reach[i]) return i; - assert(false); + vector can_reach(t+1); + can_reach[0] = true; + for (int x : w) { + for (int i = t-x; i >= 0; --i) { + if (can_reach[i]) can_reach[i+x] = true; + } + } + for (int i = t;; i--) + if (can_reach[i]) return i; + assert(false); } int main() { - const int MAX_N = 10; - const int MAX_W = 10; - const int iters = 1000000; - rep(it,0,iters) { - int n = rand() % MAX_N; - int maxw = rand() % MAX_W + 1; - vi w(n); - rep(i,0,n) - w[i] = rand()%(maxw+1); - int t = rand() % (MAX_N*maxw); - assert(naive(w,t) == knapsack(w,t)); - } - cout<<"Tests passed!"<