Project Euler

Triangles Containing the Origin II

Define: x_n = (1248n mod 32323) - 16161 y_n = (8421n mod 30103) - 15051 P_n = {(x_1, y_1), (x_2, y_2),..., (x_n, y_n)} Let C(n) be the number of triangles whose vertices are in P_n which contain th...

Source sync Apr 19, 2026

Problem #0456

Level Level 23

Solved By 487

Languages C++, Python

Answer 333333208685971546

Length 436 words

geometrymodular_arithmeticlinear_algebra

Define:

For example, .

Let be the number of triangles whose vertices are in which contain the origin in the interior.

Examples:

Find .

Problem 456: Triangles Containing the Origin II

Mathematical Analysis

Triangle Contains Origin Test

A triangle with vertices A, B, C contains the origin if and only if the origin lies on the same side of each edge. Using cross products, the origin is inside triangle (A, B, C) if the signs of:

A x B (cross product z-component: a_x * b_y - a_y * b_x)
B x C
C x A

are all the same (all positive or all negative).

Efficient Counting via Angular Sorting

Instead of checking all O(n^3) triangles, we use the following approach:

Compute the angle theta_i = atan2(y_i, x_i) for each point.
Sort points by angle.
For each point P_i, count how many points lie in the semicircle opposite to P_i (i.e., angles in [theta_i + pi, theta_i + 2*pi)).
The number of triangles NOT containing the origin can be computed from these counts, then subtracted from C(n,3).

A triangle contains the origin if and only if the origin cannot be separated from the triangle by a line through the origin. Equivalently, the three angular sectors from the origin to each pair of vertices must each be less than pi.

A triangle does NOT contain the origin iff there exists a halfplane through the origin containing all three vertices. For each point i, let f(i) be the number of points in the open semicircle starting at point i going counterclockwise for pi radians. Then:

Number of “bad” triangles = sum over i of C(f(i), 2)

And: C(n) = C(n, 3) - sum_i C(f(i), 2)

where we use a two-pointer / sliding window on the sorted angles.

Editorial

We generate all points (x_i, y_i) for i = 1 to n. We then compute angles and sort. Finally, use two-pointer technique to count f(i) for each point.

Pseudocode

Generate all points (x_i, y_i) for i = 1 to n
Compute angles and sort
Use two-pointer technique to count f(i) for each point
Compute total = C(n,3) - sum of C(f(i), 2)

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

Time: O(n log n) for sorting, O(n) for two-pointer sweep.
Space: O(n).

Answer

\boxed{333333208685971546}

C++ project_euler/problem_456/solution.cpp

#include <bits/stdc++.h>
using namespace std;

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

    auto solve = [](int N) -> long long {
        // Generate points using recurrence
        vector<long long> px(N), py(N);
        long long xv = 1248, yv = 8421;
        for(int i = 0; i < N; i++){
            px[i] = xv - 16161;
            py[i] = yv - 15051;
            xv = (xv * 1248) % 32323;
            yv = (yv * 8421) % 30103;
        }

        // Sort by angle using exact arithmetic
        // angle in [-pi, pi), map to [0, 2pi) using half-planes
        // Upper half (y > 0) or (y == 0, x > 0): first half [0, pi)
        // Lower half (y < 0) or (y == 0, x < 0): second half [pi, 2pi)
        auto half = [](long long x, long long y) -> int {
            if(y > 0) return 0;
            if(y < 0) return 1;
            if(x > 0) return 0;
            return 1; // y == 0, x < 0 (or x == 0 which shouldn't happen)
        };

        // Compare angles: return true if a comes before b in [0, 2pi) order
        auto angleLess = [&](int a, int b) -> bool {
            int ha = half(px[a], py[a]);
            int hb = half(px[b], py[b]);
            if(ha != hb) return ha < hb;
            // Same half: use cross product
            long long cross = px[a] * py[b] - py[a] * px[b];
            return cross > 0;
        };

        vector<int> idx(N);
        iota(idx.begin(), idx.end(), 0);
        sort(idx.begin(), idx.end(), [&](int a, int b){
            return angleLess(a, b);
        });

        vector<long long> sx(N), sy(N);
        for(int i = 0; i < N; i++){
            sx[i] = px[idx[i]];
            sy[i] = py[idx[i]];
        }

        // Check if point j is strictly less than pi ahead of point i (in sorted circular order)
        // Two points have angular difference exactly pi iff they are anti-parallel:
        // cross(i,j) = 0 and dot(i,j) < 0
        // Angular difference < pi iff:
        //   - half(i) == half(j) and cross(i,j) > 0 [same half, j ahead of i within half]
        //     but this is angle < pi within one half
        //   - half(i) != half(j) and... we're wrapping
        //
        // Easier: point j is within pi counterclockwise from i iff cross(i,j) > 0.
        // If cross(i,j) == 0, then either same direction (angle = 0) or anti-parallel (angle = pi).
        //   dot > 0 => same direction => angle = 0
        //   dot < 0 => anti-parallel => angle = pi
        //
        // For counting, we want f(i) = # points in half-open [angle_i, angle_i + pi).
        // "j at same angle as i" has angular difference 0, so included.
        // "j at angle_i + pi" has angular difference pi, so excluded.
        //
        // In sorted order (starting from i going forward through the circular array):
        // - Points at same angle: cross(i,j)=0, dot>0 => included
        // - Points with cross(i,j)>0 => included (angle in (0,pi))
        // - Points with cross(i,j)=0, dot<0 => excluded (angle = pi)
        // - Points with cross(i,j)<0 => excluded (angle > pi)
        //
        // BUT: in sorted circular order, we go from angle_i counterclockwise through 2pi.
        // The transition from "included" to "excluded" happens when we hit angle_i + pi.
        // In sorted order: first same-angle points, then cross>0 points, then cross=0/dot<0, then cross<0.

        long long bad = 0;
        int j = 0;
        for(int i = 0; i < N; i++){
            if(j <= i) j = i + 1;
            while(j < i + N){
                int jj = j % N;
                long long cross = sx[i] * sy[jj] - sy[i] * sx[jj];
                if(cross > 0){
                    j++;
                    continue;
                }
                if(cross == 0){
                    long long dot = sx[i] * sx[jj] + sy[i] * sy[jj];
                    if(dot > 0){
                        // same direction, angle difference = 0, include
                        j++;
                        continue;
                    }
                }
                break;
            }
            long long fi = (long long)(j - i - 1);
            bad += fi * (fi - 1) / 2;
        }

        long long total = (long long)N * (N - 1) / 2;
        total = total * (N - 2) / 3;
        return total - bad;
    };

    printf("C(8) = %lld\n", solve(8));
    printf("C(600) = %lld\n", solve(600));
    printf("C(40000) = %lld\n", solve(40000));
    printf("C(2000000) = %lld\n", solve(2000000));

    return 0;
}

Python project_euler/problem_456/solution.py

import math

def solve():
    N = 2000000

    # Generate points and compute angles
    angles = []
    for i in range(1, N + 1):
        xi = (1248 * i % 32323) - 16161
        yi = (8421 * i % 30103) - 15051
        angles.append(math.atan2(yi, xi))

    angles.sort()

    # Duplicate for circular handling
    ang = angles + [a + 2 * math.pi for a in angles]

    # Two-pointer: count f(i) = points in open semicircle
    bad = 0
    j = 0
    for i in range(N):
        if j <= i:
            j = i + 1
        while j < 2 * N and ang[j] - ang[i] < math.pi - 1e-9:
            j += 1
        fi = j - i - 1
        bad += fi * (fi - 1) // 2

    total = N * (N - 1) * (N - 2) // 6
    answer = total - bad

    print(answer)

solve()