p440.cpp
======================
#include <bits/stdc++.h>
// #include <coding_library/geometry/point.hpp>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

using coord_t = double;

struct Point {
    static constexpr coord_t eps = 1e-9;
    static inline const coord_t PI = acos((coord_t)-1.0);

    coord_t x, y;
    Point(coord_t x = 0, coord_t y = 0) : x(x), y(y) {}

    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }
    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }
    Point operator*(coord_t c) const { return Point(x * c, y * c); }
    Point operator/(coord_t c) const { return Point(x / c, y / c); }

    coord_t operator*(const Point& p) const { return x * p.x + y * p.y; }
    coord_t operator^(const Point& p) const { return x * p.y - y * p.x; }

    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; }
    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 ? x > p.x : y > p.y;
    }
    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 ? x > p.x : y >= p.y;
    }

    coord_t norm2() const { return x * x + y * y; }
    coord_t norm() const { return sqrt(norm2()); }
    coord_t angle() const { return atan2(y, x); }

    Point rotate(coord_t a) const {
        return Point(x * cos(a) - y * sin(a), x * sin(a) + y * cos(a));
    }

    Point perp() const { return Point(-y, x); }
    Point unit() const { return *this / norm(); }
    Point normal() const { return perp().unit(); }
    Point project(const Point& p) const {
        return *this * (*this * p) / norm2();
    }
    Point reflect(const Point& p) const {
        return *this * 2 * (*this * p) / norm2() - p;
    }

    friend ostream& operator<<(ostream& os, const Point& p) {
        return os << p.x << ' ' << p.y;
    }
    friend istream& operator>>(istream& is, Point& p) {
        return is >> p.x >> p.y;
    }

    friend int ccw(const Point& a, const Point& b, const Point& c) {
        coord_t v = (b - a) ^ (c - a);
        if(-eps <= v && v <= eps) {
            return 0;
        } else if(v > 0) {
            return 1;
        } else {
            return -1;
        }
    }

    friend bool point_on_segment(
        const Point& a, const Point& b, const Point& p
    ) {
        return ccw(a, b, p) == 0 && p.x >= min(a.x, b.x) - eps &&
               p.x <= max(a.x, b.x) + eps && p.y >= min(a.y, b.y) - eps &&
               p.y <= max(a.y, b.y) + eps;
    }

    friend bool point_in_triangle(
        const Point& a, const Point& b, const Point& c, const Point& p
    ) {
        int d1 = ccw(a, b, p);
        int d2 = ccw(b, c, p);
        int d3 = ccw(c, a, p);
        return (d1 >= 0 && d2 >= 0 && d3 >= 0) ||
               (d1 <= 0 && d2 <= 0 && d3 <= 0);
    }

    friend Point line_line_intersection(
        const Point& a1, const Point& b1, const Point& a2, const Point& b2
    ) {
        return a1 +
               (b1 - a1) * ((a2 - a1) ^ (b2 - a2)) / ((b1 - a1) ^ (b2 - a2));
    }

    friend bool collinear(const Point& a, const Point& b) {
        return abs(a ^ b) < eps;
    }

    friend Point circumcenter(const Point& a, const Point& b, const Point& c) {
        Point mid_ab = (a + b) / 2.0;
        Point mid_ac = (a + c) / 2.0;
        Point perp_ab = (b - a).perp();
        Point perp_ac = (c - a).perp();
        return line_line_intersection(
            mid_ab, mid_ab + perp_ab, mid_ac, mid_ac + perp_ac
        );
    }

    friend coord_t arc_area(
        const Point& center, coord_t r, const Point& p1, const Point& p2
    ) {
        coord_t theta1 = (p1 - center).angle();
        coord_t theta2 = (p2 - center).angle();
        if(theta2 < theta1 - eps) {
            theta2 += 2 * PI;
        }

        coord_t d_theta = theta2 - theta1;
        coord_t cx = center.x, cy = center.y;
        coord_t area = r * cx * (sin(theta2) - sin(theta1)) -
                       r * cy * (cos(theta2) - cos(theta1)) + r * r * d_theta;
        return area / 2.0;
    }

    friend vector<Point> intersect_circles(
        const Point& c1, coord_t r1, const Point& c2, coord_t r2
    ) {
        Point d = c2 - c1;
        coord_t dist = d.norm();

        if(dist > r1 + r2 + eps || dist < abs(r1 - r2) - eps || dist < eps) {
            return {};
        }

        coord_t a = (r1 * r1 - r2 * r2 + dist * dist) / (2 * dist);
        coord_t h_sq = r1 * r1 - a * a;
        if(h_sq < -eps) {
            return {};
        }
        if(h_sq < 0) {
            h_sq = 0;
        }
        coord_t h = sqrt(h_sq);

        Point mid = c1 + d.unit() * a;
        Point perp_dir = d.perp().unit();

        if(h < eps) {
            return {mid};
        }
        return {mid + perp_dir * h, mid - perp_dir * h};
    }

    friend optional<Point> intersect_ray_segment(
        const Point& ray_start, const Point& ray_through, const Point& seg_a,
        const Point& seg_b
    ) {
        Point ray_dir = ray_through - ray_start;
        if(ray_dir.norm2() < Point::eps) {
            return {};
        }
        Point seg_dir = seg_b - seg_a;
        coord_t denom = ray_dir ^ seg_dir;
        if(fabs(denom) < eps) {
            return {};
        }
        coord_t t = ((seg_a - ray_start) ^ seg_dir) / denom;
        if(t < eps) {
            return {};
        }
        coord_t s = ((seg_a - ray_start) ^ ray_dir) / denom;
        if(s < eps || s > 1 - eps) {
            return {};
        }
        return ray_start + ray_dir * t;
    }
};

struct Interval {
    double lo, hi;
};

vector<Interval> allowed;

void precompute_intervals() {
    static const double PI = acos(-1.0);
    static const double DEG = PI / 180.0;

    for(int k = 1; k <= 9; k++) {
        double a_lo_deg = 90.0 / (k + 1);
        double a_hi_deg = 90.0 / k;

        double cos_hi = cos(a_lo_deg * DEG);
        double cos_lo = cos(a_hi_deg * DEG);

        int m_min = (int)floor((cos_lo * 1000 - k) / 10) - 1;
        int m_max = (int)ceil((cos_hi * 1000 - k) / 10) + 1;

        for(int m = m_min; m <= m_max; m++) {
            double c_lo = max((10 * m + k) / 1000.0, cos_lo);
            double c_hi = min((10 * m + k + 1) / 1000.0, cos_hi);

            if(c_lo >= c_hi) {
                continue;
            }

            double a_int_lo = acos(min(c_hi, 1.0));
            double a_int_hi = acos(max(c_lo, 0.0));

            a_int_lo = max(a_int_lo, a_lo_deg * DEG);
            a_int_hi = min(a_int_hi, a_hi_deg * DEG);

            if(a_int_lo >= a_int_hi) {
                continue;
            }

            allowed.push_back({a_int_lo - 1e-9, a_int_hi + 1e-9});
        }
    }
}

int n;
vector<Point> pts;

void read() {
    cin >> n;
    pts.resize(n);
    for(auto& p: pts) {
        cin >> p;
    }
}

void solve() {
    // It's trivial to solve this problem in O(N^3), by trying all possible
    // positions of the three people, and then checking the angle condition.
    // However, this might end up being a bit too slow, even with optimizations.
    // Instead we can notice that the condition we want is that A = angle(QPR)
    // satisfies floor(90/A) = floor(cos(A) * 1000) % 10, and that there aren't
    // that many contiguous ranges with such angles. In particular, there are
    // only around ~100 such angle ranges. We will denote this count by K, and
    // we will precompute these.
    //
    // This immediately gives us a fairly straight forward solution that runs in
    // O(N^2 log N + N^2K). We will simply choose each possible "middle" point
    // P, and then sort by angle all other points around it. Afterwards, we can
    // make make O(K) sweeps in O(N) each that would go through the ranges we
    // defined earlier, but we should be careful about is the floating point
    // numbers, so one way to do this is to slightly relax the ranges, and then
    // trim them by explicitly checking the triplets when we do each of the
    // sweeps.

    static const double PI = Point::PI;
    static const double TWO_PI = 2.0 * PI;
    int num_intervals = (int)allowed.size();

    static const double DEG = PI / 180.0;

    double cos2_bound[11];
    for(int k = 1; k <= 10; k++) {
        double c = cos(90.0 / k * DEG);
        cos2_bound[k] = c * c;
    }

    auto check = [&](const Point& dj, int64_t L1, const Point& dk,
                     int64_t L2) -> bool {
        int64_t dot = (int64_t)(dj.x * dk.x + dj.y * dk.y);
        if(dot <= 0) {
            return false;
        }

        using i128 = __int128;
        i128 D = (i128)1000 * dot;
        i128 D2 = D * D;
        i128 L1L2 = (i128)L1 * L2;

        int F = (int)floor(1000.0 * dot / sqrt((double)L1 * L2));
        while(F > 0 && (i128)F * F * L1L2 > D2) {
            F--;
        }
        while((i128)(F + 1) * (F + 1) * L1L2 <= D2) {
            F++;
        }
        if(F >= 1000) {
            return false;
        }

        int rhs = F % 10;

        double cos2_a = (double)dot * dot / ((double)L1 * L2);
        int lhs = 0;
        for(int k = 9; k >= 1; k--) {
            if(cos2_a >= cos2_bound[k] - 1e-12) {
                if(cos2_a < cos2_bound[k + 1] + 1e-12) {
                    lhs = k;
                }
                break;
            }
        }

        return lhs > 0 && lhs == rhs;
    };

    int64_t answer = 0;
    for(int i = 0; i < n; i++) {
        int m = n - 1;
        vector<int> idx_map(m);
        vector<double> angles(m);
        vector<int64_t> norm2s(m);
        vector<Point> deltas(m);
        int idx = 0;
        for(int j = 0; j < n; j++) {
            if(j == i) {
                continue;
            }
            Point d = pts[j] - pts[i];
            idx_map[idx] = j;
            angles[idx] = d.angle();
            norm2s[idx] = (int64_t)(d.x * d.x + d.y * d.y);
            deltas[idx] = d;
            idx++;
        }

        vector<int> order(m);
        iota(order.begin(), order.end(), 0);
        sort(order.begin(), order.end(), [&](int a, int b) {
            return angles[a] < angles[b];
        });

        vector<double> ext(2 * m);
        vector<int> ext_ord(2 * m);
        for(int j = 0; j < m; j++) {
            ext[j] = angles[order[j]];
            ext[j + m] = angles[order[j]] + TWO_PI;
            ext_ord[j] = order[j];
            ext_ord[j + m] = order[j];
        }

        vector<int> ptr_lo(num_intervals, 0);
        vector<int> ptr_hi(num_intervals, 0);

        for(int j = 0; j < m; j++) {
            double base = ext[j];
            int oj = ext_ord[j];
            for(int k = 0; k < num_intervals; k++) {
                double lo = base + allowed[k].lo;
                double hi = base + allowed[k].hi;

                while(ptr_lo[k] < 2 * m && ext[ptr_lo[k]] < lo) {
                    ptr_lo[k]++;
                }
                while(ptr_hi[k] < 2 * m && ext[ptr_hi[k]] < hi) {
                    ptr_hi[k]++;
                }

                int L = ptr_lo[k];
                int R = ptr_hi[k];

                while(L < R) {
                    int op = ext_ord[L];
                    if(check(deltas[oj], norm2s[oj], deltas[op], norm2s[op])) {
                        break;
                    }
                    L++;
                }

                while(L < R) {
                    int op = ext_ord[R - 1];
                    if(check(deltas[oj], norm2s[oj], deltas[op], norm2s[op])) {
                        break;
                    }
                    R--;
                }

                if(L < R) {
                    answer += (R - L);
                }
            }
        }
    }

    cout << answer << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    precompute_intervals();

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        solve();
    }

    return 0;
}

=================
p440.in1
======================
10
628 1
17 207
176 1
16 -5
161 0
-1 56
17 83
1 5
15 1
18 101

=================
p440.ans1
======================
15

=================
statement.txt
======================
440. Moles and Holes
Time limit per test: 2.25 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard



N moles live in the cottage of Mr. Tapacob. Every mole lives in its own hole. And Mr. Tapacob is not glad to see them. So he forced them out. But there are three strongest moles, who don't want to leave. They are Shureg, Ruslan and Ruslan. And you should not distinguish Ruslans. Sometimes these three moles poke out of holes to look around. But only one of them, the strangest mole (Shureg), is not blind. And when they peek out from hole, Shureg looks at two other moles. But he is a mole and his vision is strange, therefore he can't see in any angle. The angle between two other moles (let's call it A) must be sharp, and an integer part of value 90 / A (in degrees) must be equal to the third digit after decimal point of decimal representation of .

You are given coordinates of holes. Mr. Tapacob wants to know how many ways are there for moles to peek out from holes such that Shureg can see two other moles. Any mole can get out at any of the holes, but the three holes which they get from must be all different.

Input
First line of input contains integer N (3 ≤ N ≤ 800). Next N lines contain coordinates of holes; each line consists of two integers separated by a space. Coordinates do not exceed 1000 by absolute value. No two holes coincide.

Output
Output one integer — the number of ways.

Example(s)
sample input
sample output
10
628 1
17 207
176 1
16 -5
161 0
-1 56
17 83
1 5
15 1
18 101
15

=================
