#include #include #include // <-- Begin of Algorithm Implementation/Geometry/Convex hull/Monotone chain --> // https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#C.2B.2B // Implementation of Andrew's monotone chain 2D convex hull algorithm. // Asymptotic complexity: O(n log n). // Practical performance: 0.5-1.0 seconds for n=1000000 on a 1GHz machine. #include #include using namespace std; typedef long coord_t; // coordinate type typedef long long coord2_t; // must be big enough to hold 2*max(|coordinate|)^2 struct Point { Point() : x(0), y(0), r2(0) {} Point(long long xx, long long yy) : x(xx), y(yy), r2(xx*xx + yy*yy) {} coord_t x, y; coord2_t r2; bool operator <(const Point &p) const { return x < p.x || (x == p.x && y < p.y); } bool operator ==(const Point &p) const { return x == p.x && y == p.y; } bool operator !=(const Point &p) const { return x != p.x || y != p.y; } Point operator +(const Point &p) const { return Point(this->x + p.x, this->y + p.y); } Point operator -(const Point &p) const { return Point(this->x - p.x, this->y - p.y); } }; // 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product. // Returns a positive value, if OAB makes a counter-clockwise turn, // negative for clockwise turn, and zero if the points are collinear. coord2_t cross(const Point &O, const Point &A, const Point &B) { return (A.x - O.x) * (B.y - O.y) - (A.y - O.y) * (B.x - O.x); } // Returns a list of points on the convex hull in counter-clockwise order. // Note: the last point in the returned list is the same as the first one. vector convex_hull(vector P) { int n = P.size(), k = 0; if (n == 1) return P; vector H(2*n); // Sort points lexicographically sort(P.begin(), P.end()); // Build lower hull for (int i = 0; i < n; ++i) { while (k >= 2 && cross(H[k-2], H[k-1], P[i]) <= 0) k--; H[k++] = P[i]; } // Build upper hull for (int i = n-2, t = k+1; i >= 0; i--) { while (k >= t && cross(H[k-2], H[k-1], P[i]) <= 0) k--; H[k++] = P[i]; } H.resize(k-1); return H; } // <-- End of Algorithm Implementation/Geometry/Convex hull/Monotone chain --> // scalar product coord2_t scalar(Point &a, Point &b) { return a.x * b.x + a.y * b.y; } // Compare points by norm struct QueueCompare { bool operator()(const Point& lhs, const Point& rhs) const { return lhs.r2 > rhs.r2; } }; // Walker - enumerate points x>=0 and y>=0 in increasing norm class Walker { public: Walker() { s = 1; queue.push(Point(1,0)); } Point Next() { while (queue.top().r2 >= s*s) { queue.push(Point(0, s)); s++; } Point p = queue.top(); queue.pop(); queue.push(Point(p.x+1, p.y)); return p; } private: long long s; priority_queue, QueueCompare> queue; }; vector hull; // convex hull of all previous points Point penultimate; // penultimate point Point last; // last point bool ok(const Point &p) { for (size_t h=0; h= 0) { return false; } } } return true; } Point compute(int n) { switch(n) { case 1: return Point(0,0); case 2: return Point(1,0); default: Point prev = last - penultimate; Point best; coord2_t r2 = 0; coord2_t bestScal = 0; Walker walker; for (;;) { Point delta = walker.Next(); if (r2 && r2 < delta.r2) { break; } for (int r=0; r<4; r++) { coord2_t scal = scalar(prev, delta); if ( r2==0 || bestScal < scal ) { Point o = last + delta; if (cross(last, penultimate, o) < 0 && ok(o) ) { r2 = delta.r2; best = o; bestScal = scal; } } // rotate by 90 degrees delta = Point(-delta.y, +delta.x); } } return best; } } int main() { printf("background red\n"); for (int n=1; n<=500; n++) { Point p = compute(n); hull.push_back(p); hull = convex_hull(hull); printf("%d %d\n", n, p.y); fflush(stdout); penultimate = last; last = p; } return 0; }