## 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. Detailed editorial

### 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.

The condition is:

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

Since `0 < A < 90`, we have:

```text
floor(90 / A) >= 1
```

The right-hand side is a decimal digit, so it is between `0` and `9`.

Therefore:

```text
floor(90 / A) in {1, 2, ..., 9}
```

So the possible value of `floor(90 / A)` is small.

---

### Observation 2: precompute valid angle intervals

Let:

```text
k = floor(90 / A)
```

Then:

```text
k <= 90 / A < k + 1
```

which implies:

```text
90 / (k + 1) < A <= 90 / k
```

for `k = 1..9`.

We also need:

```text
floor(1000 * cos(A)) mod 10 = k
```

That means there exists some integer `m` such that:

```text
10m + k <= floor(1000 * cos(A)) < 10m + k + 1
```

Equivalently:

```text
(10m + k) / 1000 <= cos(A) < (10m + k + 1) / 1000
```

Since `cos(A)` is monotonic decreasing on `(0, 90)`, each such cosine interval corresponds to one angle interval.

For every `k = 1..9`, we intersect:

```text
90 / (k + 1) < A <= 90 / k
```

with all cosine-digit intervals where the third decimal digit is `k`.

This produces only about `100` valid angle intervals.

Call the number of intervals `K`.

---

### Observation 3: process each Shureg hole independently

Fix a hole `P` as Shureg's hole.

For every other point `X`, compute the polar angle of vector:

```text
X - P
```

Sort all other points by this polar angle.

Now, for a fixed first ray `P -> Q`, we need to find all other rays `P -> R` such that the counter-clockwise angular difference from `Q` to `R` lies in one of the valid intervals.

Because the valid angle is acute, i.e. `< 90`, for every unordered pair `{Q, R}` exactly one of the two directed angular differences is the actual small angle. Therefore counting directed differences in valid intervals counts each unordered Ruslan pair exactly once.

---

### Observation 4: duplicate the angle array

Angles wrap around at `2*pi`.

To handle this easily, duplicate the sorted angle array:

```text
angles[i + m] = angles[i] + 2*pi
```

where `m = N - 1`.

Now for every base angle `angles[j]`, all possible forward angular differences can be searched in the interval:

```text
[j + 1, j + m)
```

inside the duplicated array.

---

### Observation 5: two-pointers for every valid interval

For each valid angle interval `[lo, hi]`, and for each base angle `base`, we need all points with angle in:

```text
[base + lo, base + hi)
```

Since base angles are processed in increasing order, the interval endpoints also move monotonically.

Therefore, for each valid interval, maintain two pointers:

```text
ptr_lo
ptr_hi
```

These pointers move only forward, so processing all bases for one interval costs `O(N)`.

For one fixed Shureg hole, this gives:

```text
O(N log N + N * K)
```

because of sorting plus sweeping.

For all holes:

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

With `N <= 800` and `K ~= 100`, this is fast enough.

---

### Exact checking

Floating-point intervals are slightly enlarged to avoid missing answers due to precision errors.

However, when actually counting candidate pairs, the solution verifies the condition exactly enough using integer arithmetic:

- Dot product is integer.
- Squared lengths are integers.
- `floor(1000 * cos(A))` is computed with correction loops using integer comparisons.
- `floor(90 / A)` is determined using precomputed cosine thresholds.

This prevents incorrect counting near interval boundaries.

---

### Complexity

Let `K` be the number of valid angle intervals, about `100`.

For each center:

```text
O(N log N + N*K)
```

For all centers:

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

Memory usage:

```text
O(N + K)
```

---

## 3. C++ Solution

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

---

## 4. Python solution

```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()
```

---

## 5. Compressed editorial

The validity of a triple depends only on the angle `A` at Shureg's hole.

Since:

```text
floor(90/A)
```

must equal a decimal digit, its value can only be `1..9`.

For each `k = 1..9`, the condition:

```text
floor(90/A) = k
```

gives:

```text
90/(k+1) < A <= 90/k
```

The digit condition:

```text
floor(1000*cos(A)) % 10 = k
```

splits this range into small cosine intervals. Convert them back to angle intervals and precompute all valid intervals. There are only about `100`.

Then fix every hole as Shureg's hole. Sort all other holes by polar angle around it. Duplicate the angle array with `+2*pi` to handle wrap-around.

For each base ray and for every valid interval `[lo, hi]`, use two pointers to find all rays whose angular difference lies in:

```text
[base + lo, base + hi)
```

Because the target angle is acute, every unordered pair of Ruslan holes is counted exactly once.

Use an exact final check based on dot products and squared lengths to avoid floating-point boundary errors.

Complexity:

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

where `K ~= 100`, fast enough for `N <= 800`.
