<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

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 (with A measured in degrees) must be equal to the third digit after the decimal point of the decimal representation of cos A.

Here A is the angle of the triangle at Shureg's hole (the vertex from which he looks), i.e. the angle between the directions to the two other moles. The condition is that A is acute (0 < A < 90 degrees) and that floor(90 / A) equals the third digit after the decimal point of cos A; that is, floor(90 / A) = floor(1000 * cos A) mod 10. Note that since 0 < A < 90, the value floor(90 / A) is at least 1, so triplets with the third decimal digit of cos A equal to 0 never count.

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

<|response|>
## 1. Abridged problem statement

Given `N` distinct holes with integer coordinates (`3 ≤ N ≤ 800`), count the number of ways to choose three distinct holes for three moles:

- Shureg is distinguishable and stands at one chosen hole.
- The two Ruslans are indistinguishable and stand at the other two holes.

Let `A` be the angle at Shureg's hole between the directions to the two Ruslans, measured in degrees.

The choice is valid if:

```text
0 < A < 90
floor(90 / A) = floor(1000 * cos(A)) mod 10
```

Output the number of valid choices.

Because the Ruslans are indistinguishable, for each Shureg hole we count unordered pairs of other holes.

---

## 2. Key observations

### Observation 1: the angle condition depends only on `A`

For a fixed Shureg position `P`, and two other holes `Q` and `R`, only the angle

```text
A = angle(QPR)
```

matters.

Since `0 < A < 90`, we have `floor(90 / A) >= 1`. The right-hand side is a decimal digit `0..9`. Therefore `floor(90 / A)` can only be `1..9`.

---

### Observation 2: precompute valid angle intervals

Let `k = floor(90 / A)`. Then `90 / (k + 1) < A <= 90 / k`.

We also need `floor(1000 * cos(A)) mod 10 = k`, which means `cos(A)` must fall in one of several small intervals. Since `cos(A)` is monotone decreasing on `(0, 90)`, each cosine interval maps to an angle interval.

For `k = 1..9`, intersecting the angle range for that `k` with all cosine-digit intervals where the last digit equals `k` yields about `100` valid angle intervals total.

---

### Observation 3: process each Shureg hole independently with sorted angles and two pointers

Fix a hole `P` as Shureg's hole. Compute the polar angle of every other hole relative to `P`, then sort. For each valid angle interval `[lo, hi]`, use two pointers scanning the sorted (duplicated) angle array to count all second rays `Q -> R` whose angular difference from a first ray `P -> Q` lies in `[lo, hi)`. Because the angle is acute, each unordered Ruslan pair is counted exactly once.

Slightly enlarge the floating-point intervals and use an exact integer-arithmetic check near the boundaries to avoid precision errors.

Complexity:

```text
O(N^2 log N + N^2 K)
```

where `K ~= 100`.

---

## 3. Full solution approach

1. Call `precompute_intervals()` (called from main) to build all valid angle intervals.
2. For each hole `i` as Shureg:
   - Compute vector `d = pts[j] - pts[i]`, polar angle, squared length, for all `j != i`.
   - Sort by polar angle.
   - Duplicate the angle array with `+2*pi` for wrap-around.
   - For each base ray `j` and each valid interval, advance two pointers to find the candidate range.
   - Trim false positives at boundaries using exact integer arithmetic.
   - Add the count to the answer.
3. Output the answer.

---

## 4. C++ implementation

```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;
}
```

---

## 5. Python implementation

```python
import sys
import math


# ------------------------------------------------------------
# Precompute valid angle intervals
# ------------------------------------------------------------

def precompute_intervals():
    # List of valid angle intervals in radians.
    allowed = []

    # Degree to radian conversion.
    DEG = math.pi / 180.0

    # k = floor(90 / A), and k must be a digit 1..9.
    for k in range(1, 10):
        # From floor(90/A) = k:
        # 90/(k+1) < A <= 90/k.
        a_lo_deg = 90.0 / (k + 1)
        a_hi_deg = 90.0 / k

        # Cosine bounds for this angle range.
        # cos is decreasing on [0, 90 degrees].
        cos_hi = math.cos(a_lo_deg * DEG)
        cos_lo = math.cos(a_hi_deg * DEG)

        # Need floor(1000*cos(A)) % 10 == k.
        #
        # Therefore:
        # 1000*cos(A) in [10m+k, 10m+k+1)
        m_min = int(math.floor((cos_lo * 1000.0 - k) / 10.0)) - 1
        m_max = int(math.ceil((cos_hi * 1000.0 - k) / 10.0)) + 1

        # Try every relevant block.
        for m in range(m_min, m_max + 1):
            # Intersect cosine interval with this decimal-digit block.
            c_lo = max((10 * m + k) / 1000.0, cos_lo)
            c_hi = min((10 * m + k + 1) / 1000.0, cos_hi)

            # Empty intersection.
            if c_lo >= c_hi:
                continue

            # Convert cosine interval back to angle interval.
            # Larger cosine means smaller angle.
            a_int_lo = math.acos(min(c_hi, 1.0))
            a_int_hi = math.acos(max(c_lo, 0.0))

            # Intersect with k's angle range.
            a_int_lo = max(a_int_lo, a_lo_deg * DEG)
            a_int_hi = min(a_int_hi, a_hi_deg * DEG)

            # Empty interval.
            if a_int_lo >= a_int_hi:
                continue

            # Slightly enlarge interval to avoid floating-point misses.
            allowed.append((a_int_lo - 1e-9, a_int_hi + 1e-9))

    return allowed


# ------------------------------------------------------------
# Main solver
# ------------------------------------------------------------

def solve():
    data = sys.stdin.buffer.read().split()

    n = int(data[0])

    pts = []

    pos = 1
    for _ in range(n):
        x = int(data[pos])
        y = int(data[pos + 1])
        pos += 2
        pts.append((x, y))

    # Precomputed allowed angular intervals.
    allowed = precompute_intervals()
    num_intervals = len(allowed)

    PI = math.pi
    TWO_PI = 2.0 * PI
    DEG = PI / 180.0

    # cos2_bound[k] = cos(90/k degrees)^2.
    cos2_bound = [0.0] * 11
    for k in range(1, 11):
        c = math.cos(90.0 / k * DEG)
        cos2_bound[k] = c * c

    # --------------------------------------------------------
    # Exact-ish validity check for two vectors.
    # --------------------------------------------------------
    def check(dx1, dy1, len1, dx2, dy2, len2):
        # Dot product.
        dot = dx1 * dx2 + dy1 * dy2

        # Non-acute angle.
        if dot <= 0:
            return False

        # We need floor(1000*cos(A)).
        #
        # cos(A) = dot / sqrt(len1 * len2)
        prod_len = len1 * len2
        D = 1000 * dot
        D2 = D * D

        # Initial floating-point estimate.
        F = int(math.floor(1000.0 * dot / math.sqrt(prod_len)))

        # Correct possible overestimate.
        while F > 0 and F * F * prod_len > D2:
            F -= 1

        # Correct possible underestimate.
        while (F + 1) * (F + 1) * prod_len <= D2:
            F += 1

        # F == 1000 means angle is zero, not allowed.
        if F >= 1000:
            return False

        # Third digit after the decimal point of cos(A).
        rhs = F % 10

        # Compute cos(A)^2.
        cos2_a = (dot * dot) / prod_len

        # Determine lhs = floor(90 / A) using cosine thresholds.
        lhs = 0

        for k in range(9, 0, -1):
            if cos2_a >= cos2_bound[k] - 1e-12:
                if cos2_a < cos2_bound[k + 1] + 1e-12:
                    lhs = k
                break

        return lhs > 0 and lhs == rhs

    answer = 0

    # --------------------------------------------------------
    # Try each point as Shureg's hole.
    # --------------------------------------------------------
    for center in range(n):
        cx, cy = pts[center]

        vectors = []

        # Build vectors from center to all other points.
        for j in range(n):
            if j == center:
                continue

            x, y = pts[j]

            dx = x - cx
            dy = y - cy

            angle = math.atan2(dy, dx)
            length2 = dx * dx + dy * dy

            vectors.append((angle, dx, dy, length2))

        # Sort vectors by polar angle.
        vectors.sort(key=lambda item: item[0])

        m = n - 1

        # Separate sorted arrays for speed/readability.
        angles = [v[0] for v in vectors]
        dxs = [v[1] for v in vectors]
        dys = [v[2] for v in vectors]
        lens = [v[3] for v in vectors]

        # Duplicate angles to handle circular wrap-around.
        ext_angles = angles + [a + TWO_PI for a in angles]

        # ext_idx maps duplicated positions back to original sorted vector index.
        ext_idx = list(range(m)) + list(range(m))

        # Two pointers for every interval.
        ptr_lo = [0] * num_intervals
        ptr_hi = [0] * num_intervals

        # Use each vector as the first ray.
        for j in range(m):
            base = ext_angles[j]

            dx1 = dxs[j]
            dy1 = dys[j]
            len1 = lens[j]

            # Sweep all allowed intervals.
            for interval_id, (lo_add, hi_add) in enumerate(allowed):
                lo = base + lo_add
                hi = base + hi_add

                # Advance lower pointer.
                p = ptr_lo[interval_id]
                while p < 2 * m and ext_angles[p] < lo:
                    p += 1
                ptr_lo[interval_id] = p

                # Advance upper pointer.
                q = ptr_hi[interval_id]
                while q < 2 * m and ext_angles[q] < hi:
                    q += 1
                ptr_hi[interval_id] = q

                L = p
                R = q

                # Trim invalid candidates from left boundary.
                while L < R:
                    idx2 = ext_idx[L]

                    if check(
                        dx1, dy1, len1,
                        dxs[idx2], dys[idx2], lens[idx2]
                    ):
                        break

                    L += 1

                # Trim invalid candidates from right boundary.
                while L < R:
                    idx2 = ext_idx[R - 1]

                    if check(
                        dx1, dy1, len1,
                        dxs[idx2], dys[idx2], lens[idx2]
                    ):
                        break

                    R -= 1

                # Remaining candidates are valid.
                if L < R:
                    answer += R - L

    print(answer)


if __name__ == "__main__":
    solve()
```
