Problem 165: Intersections

Project Euler Problem 165 Solution

EulerSolve provides an optimized solution for Project Euler Problem 165, Intersections, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 165 asks for the number of distinct true intersection points among 5000 line segments generated by a deterministic pseudo-random process. Each segment is determined by four consecutive values from the sequence, so the geometric input is fixed even though it looks random. A true intersection means a single point that lies strictly inside both segments. Endpoint touches are excluded, and collinear overlaps are excluded as well. The two real difficulties are deciding that strict condition correctly for every pair of segments and making sure that the same geometric point is not counted twice when different segment pairs meet there. Mathematical Approach The solution has four mathematical ingredients: generating the segments, testing whether two segments cross in their interiors, computing the crossing point exactly, and storing those points in a canonical form so duplicates collapse automatically. The Deterministic Segment Family The sequence used by the problem is $$s_0=290797,\qquad s_{n+1}=s_n^2 \bmod 50515093,$$ and then $$t_n=s_n \bmod 500.$$ The $n$-th segment is therefore $$L_n=\bigl((t_{4n-3},t_{4n-2}),(t_{4n-1},t_{4n})\bigr),\qquad 1\le n\le 5000.$$ All coordinates are integers between 0 and 499, so every segment endpoint lies on the integer lattice inside a $500\times500$ box....

Detailed mathematical approach

Problem Summary

A segment is fixed by its two endpoints, and two segments in the plane have no common point, exactly one, or infinitely many. A common point $T$ of two segments is a true intersection point if it is their only common point and an interior point of both, that is, an endpoint of neither. In the example of the statement, $L_1$ runs from $(27,44)$ to $(12,32)$, $L_2$ from $(46,53)$ to $(17,62)$ and $L_3$ from $(46,70)$ to $(22,40)$: $L_2$ and $L_3$ have a true intersection point, the endpoint $(22,40)$ of $L_3$ lies on $L_1$, which does not count, and $L_1$ and $L_2$ have no common point, so there is one true intersection point. The real input has $5000$ segments made from the numbers

$$\begin{aligned}s_0&=290797,\\ s_{n+1}&=s_n^2 \bmod 50515093,\\ t_n&=s_n \bmod 500;\end{aligned}$$

segment $k$ runs from $(t_{4k-3},t_{4k-2})$ to $(t_{4k-1},t_{4k})$ for $k=1,\dots,5000$. The first four numbers are $27$, $144$, $12$ and $232$, so the first segment runs from $(27,144)$ to $(12,232)$. We must count the distinct points that are a true intersection point of at least one pair of the $5000$ segments. ^[1]

The whole solution in one sentence: for each of the $12\,497\,500$ pairs of segments, the signs of four signed areas decide whether the two segments cross at a point inside both; each crossing point is written exactly as two fractions in lowest terms, and a set of these points counts a point where several segments cross only once.

Mathematical Approach

For every pair of segments we must decide whether they cross, find where, and recognise a point found before. Integer arithmetic does all three exactly.

1. The input segments

The generator gives $s_1=629527$, $s_2=13339144$, $s_3=15552512$ and $s_4=17939732$, with the remainders $27,144,12,232$ modulo $500$. The next numbers $t_5,\dots,t_{12}$ are $46,53,117,462,246,70,322,40$, so the second and third segments run from $(46,53)$ to $(117,462)$ and from $(246,70)$ to $(322,40)$. The example segments $L_1$, $L_2$, $L_3$ are these first three segments with every coordinate taken modulo $100$. All coordinates are integers from $0$ to $499$; no segment has length $0$, and none occurs twice. The modulus $50515093=5807\cdot8699$ is a product of two primes of the form $4j+3$, as the Blum Blum Shub generator prescribes. ^[2]

2. Which side of a line a point lies on

For points $A$, $B$, $C$ put

$$\begin{aligned}\Delta(A,B,C)&=(B_x-A_x)(C_y-A_y)\\&-(B_y-A_y)(C_x-A_x),\end{aligned}$$

the cross product of $B-A$ and $C-A$, which is twice the signed area of the triangle $ABC$. It is positive when $C$ lies to the left of the line from $A$ towards $B$, negative when $C$ lies to the right, and $0$ exactly when the three points are collinear. ^[3] For $A=(0,0)$ and $B=(4,0)$, the points $(1,2)$, $(3,-1)$ and $(7,0)$ give $8$, $-4$ and $0$.

Fix $A\ne B$ and write $f(P)=\Delta(A,B,P)$. When $P_\lambda=P_0+\lambda(P_1-P_0)$ runs along a segment, $0\le\lambda\le1$,

$$f(P_\lambda)=(1-\lambda)\,f(P_0)+\lambda\,f(P_1),$$

because $f(P)$ is the cross product of the fixed vector $B-A$ with $P-A$, and a cross product is linear in each factor. So $f$ changes evenly along a segment, and $f(P_1)-f(P_0)$, the cross product of $B-A$ and $P_1-P_0$, is $0$ exactly when the two directions are parallel.

3. When two segments cross

Let $AB$ and $CD$ be segments of positive length, and put

$$\alpha=\Delta(A,B,C),\quad \beta=\Delta(A,B,D),$$

$$\gamma=\Delta(C,D,A),\quad \delta=\Delta(C,D,B).$$

Claim. The segments have a true intersection point if and only if $\alpha\beta\lt0$ and $\gamma\delta\lt0$: $C$ and $D$ lie strictly on opposite sides of the line $AB$, and $A$ and $B$ strictly on opposite sides of the line $CD$. The usual test for whether two closed segments meet uses $\le$ in place of $\lt$ and treats collinear segments separately. ^[3]

Proof. If $\alpha\beta\lt0$, then $f(P)=\Delta(A,B,P)$ changes evenly from $\alpha$ to $\beta$ along $CD$ and vanishes at exactly one point of it,

$$X=C+u\,(D-C),\quad u=\frac{\alpha}{\alpha-\beta},$$

with $0\lt u\lt1$. So $X$ is an interior point of $CD$ on the line $AB$. As $\alpha\ne\beta$, the directions are not parallel, so the two lines meet only in $X$. Likewise $\gamma\delta\lt0$ gives an interior point of $AB$ on the line $CD$, which lies on both lines and is therefore $X$ again. Thus $X$ is interior to both segments and, as the lines meet only once, it is their only common point.

Conversely, let $T$ be a true intersection point. Parallel directions would put both segments on one line through $T$, and as $T$ is interior to both, they would share a whole piece of it, infinitely many points. So $\alpha\ne\beta$. Writing $T=C+u(D-C)$ with $0\lt u\lt1$ gives $0=f(T)=(1-u)\alpha+u\beta$. Here $\alpha=0$ would force $\beta=0=\alpha$; hence $\alpha\ne0$, and $\beta=-(1-u)\alpha/u$ has the opposite sign. Exchanging the two segments gives $\gamma\delta\lt0$.

The strict signs settle every special case: parallel segments give $\alpha=\beta$, collinear ones four zeros, and an endpoint on the other line a zero. All of these occur in the input: $462$ pairs are parallel, one of them a pair of overlapping collinear segments, and $759$ other pairs meet in a single point that is an endpoint of at least one of them. A test with $\le$ would accept those $759$ pairs and the collinear pair.

4. The three segments of the statement

For $L_1=AB$ with $A=(27,44)$, $B=(12,32)$ and $L_2=CD$ with $C=(46,53)$, $D=(17,62)$ we find $\alpha=93$ and $\beta=-390$, so the line of $L_1$ separates the endpoints of $L_2$; but $\gamma=432$ and $\delta=915$ have the same sign, so both endpoints of $L_1$ lie on one side of the line of $L_2$, and there is no crossing. For $L_1$ and $L_3$ the four values are $-162$, $0$, $54$ and $-108$. The zero says that the endpoint $(22,40)$ of $L_3$ lies on the line of $L_1$, in fact on $L_1$ itself, since $(27,44)+\tfrac13(-15,-12)=(22,40)$; the pair is rejected, as the definition requires. For $L_2$ and $L_3$, with $L_2$ now as $AB$, the values are $\alpha=-493$, $\beta=593$, $\gamma=408$ and $\delta=-678$. Both products are negative, so the three segments have exactly one true intersection point.

5. The crossing point as exact fractions

The proof also locates the point: $X=C+u(D-C)$ with $u=\alpha/(\alpha-\beta)$, and likewise $X=A+t(B-A)$ with $t=\gamma/(\gamma-\delta)$; these are the parameters $t$ and $u$ of the usual formula for the intersection of two lines. ^[4] All the numbers involved are integers, so the coordinates of $X$ are fractions. For $L_2$ and $L_3$, $t=408/1086=68/181$ and $u=493/1086$, and

$$\begin{aligned}X&=(46,53)+\tfrac{68}{181}\,(-29,9)\\&=\Bigl(\tfrac{6354}{181},\tfrac{10205}{181}\Bigr),\end{aligned}$$

about $(35.105,\,56.381)$. Up to sign, both denominators $\gamma-\delta$ and $\alpha-\beta$ equal the cross product of the directions $B-A$ and $D-C$, here $1086$.

The C++ and Python programs use Cramer's rule instead. ^[5] The line through $(x_1,y_1)$ and $(x_2,y_2)$ is $a_1x+b_1y=c_1$ with $a_1=y_2-y_1$, $b_1=x_1-x_2$ and $c_1=a_1x_1+b_1y_1=x_1y_2-x_2y_1$. With the second line $a_2x+b_2y=c_2$, subtracting $b_1$ times the second equation from $b_2$ times the first removes $y$, and subtracting $a_2$ times the first from $a_1$ times the second removes $x$:

$$\begin{aligned}x&=\frac{c_1b_2-c_2b_1}{a_1b_2-a_2b_1},\\ y&=\frac{a_1c_2-a_2c_1}{a_1b_2-a_2b_1},\end{aligned}$$

and the denominator is again the cross product of the directions. For $L_2$ and $L_3$ the lines are $9x+29y=1951$ and $-30x+24y=300$, the denominator is $9\cdot24+30\cdot29=1086$, and the numerators are $38124$ and $61230$.

6. Why floating-point numbers are not enough

Counting floating-point coordinates is tempting but unsafe. Among the first $914$ segments every crossing point is found by one pair only. Segment $915$, from $(484,409)$ to $(194,110)$, then passes through the crossing point $(85826/269,\,64274/269)$ of segment $529$, from $(373,316)$ to $(163,16)$, and segment $653$, from $(366,224)$ to $(36,329)$. In double precision, computing $t$ and then $x_1+t\,(x_2-x_1)$ gives 319.0557620817844 along segment $529$ or $653$, but 319.05576208178434 along segment $915$, and the $y$-coordinates differ as well: a program that computes the pair of segments $529$ and $915$ along segment $915$ counts this one point twice. Computing every pair along its first segment happens to keep this point together, but not others: for the first $3000$ segments the resulting set of pairs of doubles holds $1\,027\,124$ entries instead of $1\,027\,115$.

Rounding to six decimals fails as well, because distinct crossing points can be extremely close: segment $349$ is crossed by segments $2351$ and $3447$ at two points less than $10^{-7}$ apart, which get the same rounded coordinates. Among the first $3000$ segments two pairs of points less than $10^{-6}$ apart merge in this way, leaving $1\,027\,113$ points. Floating point can be made safe only with an argument: in Cramer's rule the integer numerators and the denominator are exact and a division is correctly rounded, so equal fractions give equal doubles. ^[6] Different fractions $a/b\ne c/d$ with denominators up to $498002$ (Section 9) differ by $|ad-bc|/(bd)$, which is at least $1/498002^2\approx4\cdot10^{-12}$, far more than the largest gap $2^{-44}\approx6\cdot10^{-14}$ between neighbouring doubles below $512$, so they cannot round to the same double. Exact fractions need no such argument.

7. One key for every point

A fraction has many names, such as $38124/1086=6354/181$, so the programs bring each coordinate into one standard form: a positive denominator, by changing both signs if necessary, and numerator and denominator divided by their greatest common divisor. Every rational number has exactly one such form. ^[7] Indeed, let $p/q=p'/q'$ with $q,q'\gt0$, $\gcd(p,q)=1$ and $\gcd(p',q')=1$. Then $pq'=p'q$, so $q$ divides $pq'$ and, sharing no prime factor with $p$, divides $q'$; in the same way $q'$ divides $q$, so $q=q'$ and $p=p'$. Two pairs of segments therefore give the same four integers, numerator and denominator of $x$ and of $y$, exactly when they give the same point. The coordinates are reduced separately, so their denominators may differ; at the point of Section 6 both are $269$.

8. Counting points, not pairs

A point can be found by several pairs in two ways. If $m$ segments cross pairwise at one point, all $\binom m2$ pairs find it, like the three pairs among segments $529$, $653$ and $915$. And two overlapping collinear segments have no true intersection point with each other, but a third segment crossing their common piece crosses both at the same point. The input has such a pair: segment $1264$, from $(453,51)$ to $(485,83)$, and segment $2889$, from $(460,58)$ to $(402,0)$, lie on the line $y=x-402$ and overlap for $453\le x\le460$. A point also counts if it is an endpoint of a third segment, which happens at $16$ points of the input. For the first $n$ segments:

$n=1000$: $113\,853$ crossing pairs, $113\,849$ points, two of them crossings of three segments;
$n=2000$: $463\,882$ crossing pairs, $463\,866$ points, eight of them crossings of three segments;
$n=3000$: $1\,027\,166$ crossing pairs and $1\,027\,115$ points; $19$ crossings of three segments are found three times and $13$ points on the overlap of segments $1264$ and $2889$ twice, so $1\,027\,166-2\cdot19-13=1\,027\,115$.

In the whole input at most three segments cross at one point. A set of exact points absorbs all these repetitions without special treatment.

9. How large the numbers get

Coordinates lie between $0$ and $499$, so each component of a direction is at most $499$ in absolute value. Hence the four signed areas and the common denominator are at most $2\cdot499^2=498002$ in absolute value, $|c_1|=|x_1y_2-x_2y_1|$ is at most $499^2=249001$, and each numerator of Cramer's rule is at most $2\cdot249001\cdot499=248502998$ in absolute value, below $2^{31}-1=2147483647$. Even 32-bit integers would do; the programs use 64-bit integers in C++ and Java and unlimited ones in Python. The generator does need 64 bits: $s_n^2$ is less than $50515093^2\approx2.6\cdot10^{15}$, far beyond the 32-bit range but far below $2^{63}$.

10. How many pairs cross

About $23\%$ of all pairs cross: $22.8\%$ of the pairs among the first $1000$ segments, $23.2\%$ among the first $2000$ and $22.8\%$ among the first $3000$. A model with random segments explains this value, although it proves nothing about this particular input. Four points chosen independently and uniformly in a square are the corners of a convex quadrilateral with probability $25/36$, the answer to Sylvester's four-point problem for the square; ^[8] otherwise one of them lies inside the triangle of the other three, and no two segments with these endpoints cross. In the convex case two segments cross exactly when they are the diagonals, that is, when the corner opposite one endpoint of the first segment is its other endpoint; by symmetry each of the three other points is equally likely to be that corner. So two random segments cross with probability $\frac13\cdot\frac{25}{36}=\frac{25}{108}\approx23.1\%$. Sweep-line algorithms such as the Bentley–Ottmann algorithm report the $k$ crossings of $n$ segments in time $O((n+k)\log n)$ and win when $k$ is small, ^[9] but here $k$ is about $23\%$ of all pairs, so their bound becomes $O(n^2\log n)$, no better than the much simpler double loop.

How the Code Works

The C++ program ^[10] starts in main, which recognises --segments=N (default $5000$) and --skip-checkpoints and returns $1$ for an unknown option or a count below $2$; its digit parser does not guard against counts too large for an int. Unless told to skip them, run_checkpoints checks that the first generated segment is $(27,144)$–$(12,232)$ and that the three segments of the statement give exactly one true intersection point; on a failure the program returns $2$. Then solve calls generate_segments, which computes $s_0,\dots,s_{4n}$ in 64-bit arithmetic and takes the coordinates of segment $i$, counted from $0$, from $s_{4i+1},\dots,s_{4i+4}$ modulo $500$ (Section 1).

count_true_intersections runs over all pairs $i\lt j$, computes $\alpha,\beta,\gamma,\delta$ of Section 3 with cross, the function $\Delta$, as c1 to c4, and skips the pair when sgn(c1) * sgn(c2) >= 0 or sgn(c3) * sgn(c4) >= 0. Otherwise intersection_point evaluates Cramer's rule (Section 5), makes the denominator positive by changing its sign and the signs of both numerators if necessary, and divides each coordinate's numerator and denominator by their std::gcd (Section 7). The resulting RationalPoint of four integers has an operator< that compares them in order with std::tie, so a std::set, an ordered container with logarithmic insertion, keeps one copy of each point (Section 8); its size is printed.

The Python program does the same inside one function solve with nested helpers cross, sgn and intersection_point, storing points as tuples in a built-in set; it has no options and no checkpoints.

The Java program keeps the $s_n$ in an int array, squaring them as long, and tests the parameters of Section 5 instead of signs. It computes the cross product denom of the directions (dx1, dy1) and (dx2, dy2), skips parallel pairs with denom == 0, and computes tn, equal to $\gamma$, and un, equal to $-\alpha$, so that $t$ is tn/denom and $u$ is un/denom. It keeps the pair only if $0\lt t\lt1$ and $0\lt u\lt1$, checked without division in two cases for the sign of denom; this is equivalent to the test of Section 3. The coordinate $x_1+t\,(x_2-x_1)$ becomes the fraction (x1*denom + tn*dx1)/denom, $y$ likewise; each is reduced with Euclid's algorithm in gcd and given a positive denominator. Then hash4 packs the point $(a/b,\,c/d)$ into one long, with $m=1000003$, as

$$((a\cdot m+b)\cdot m+c)\cdot m+d,$$

where overflowing products wrap around, and a HashSet<Long> collects the keys. This packing is not one-to-one: $(1/3,\,1/2011)$ and $(1/2,\,1000004/2011)$ get the same key because $1000004=m+1$, so in general two points could be counted once. The sets of C++ and Python store the four integers themselves and compare them in full, so this cannot happen there. We checked that all true intersection points of this input have different keys, and the Java program prints the same count as the others.

Complexity and Verification

For $n$ segments the programs test $n(n-1)/2$ pairs, $12\,497\,500$ for $n=5000$, each with a fixed number of multiplications: eight for the four signed areas in C++ and Python, at most six for denom, tn and un in Java. Each of the $k$ crossing pairs adds two gcd computations on numbers below $2^{31}$ and one insertion into the set, which costs $O(\log k)$ comparisons in the ordered set of the C++ program and expected constant time in the hash sets of Python and Java. As $k\le n(n-1)/2$, the time is $O(n^2\log n)$ for C++ and expected $O(n^2)$ for the others, and the memory is $O(n+k)$, dominated by the set; for the first $3000$ segments it holds $1\,027\,115$ points. On a desktop computer the full run takes about $2$ seconds in C++, under $1$ second in Java and about $9$ seconds in Python.

The C++ program checks the first segment and the example of the statement before it starts. For this article all three programs were run, and they print the same count. Two enumerations written separately confirm it: a brute force straight from the definition, which solves the two parametric equations of each pair with exact fractions and treats parallel and collinear pairs on their own, and an enumeration with $X=C+u(D-C)$ from Section 3. Both agree with the programs for all $5000$ segments and with the C++ option --segments=N for $N=3$, $10$, $20$, $100$, $1000$, $2000$ and $3000$, where the first four counts are $0$, $5$, $32$ and $1112$ and the others are those of Section 8. The second enumeration also produced the special points and counts of Sections 3, 6 and 8 and showed that the Java keys do not collide. A simulation with $2\cdot10^7$ pairs of random segments in a square gives a crossing frequency of $0.2314$, close to $25/108\approx0.2315$. None of these checks is part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The crossing test with its proof, the exact crossing points, the counts for fewer segments and the floating-point experiments are derived or computed in this article.

Project Euler 165 — Intersections. The official statement defines a true intersection point as the only common point of two segments, interior to both, gives three example segments with one such point and the numbers sₙ and tₙ, whose first four tₙ are 27, 144, 12, 232, and asks how many distinct true intersection points the 5000 segments have.
Wikipedia — Blum Blum Shub. The generator of L. Blum, M. Blum and M. Shub (1986) computes xₙ₊₁ = xₙ² mod M, where M is a product of two large primes congruent to 3 modulo 4, and outputs the parity or a few low bits of each value.
cp-algorithms — Oriented area of a triangle and Check if two segments intersect. Twice the signed area is a 2 × 2 determinant, negative for a clockwise turn, positive for an anticlockwise one and zero for collinear points. Two segments that are not collinear meet when the endpoints of neither lie on the same side of the line of the other; touching counts there.
Wikipedia — Line–line intersection. The intersection point is a quotient of determinants; with parameters t and u it lies on both segments when 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, and the denominator is zero when the lines are parallel or coincident.
Wikipedia — Cramer's rule. If a₁b₂ − b₁a₂ ≠ 0, the system a₁x + b₁y = c₁, a₂x + b₂y = c₂ has the unique solution x = (c₁b₂ − b₁c₂)/(a₁b₂ − b₁a₂), y = (a₁c₂ − c₁a₂)/(a₁b₂ − b₁a₂).
Wikipedia — IEEE 754. Unless specified otherwise, an operation returns its exact result rounded by the active rounding rule, by default to nearest with ties to even.
Wikipedia — Irreducible fraction. Every rational number has a unique representation as an irreducible fraction with a positive denominator, by unique prime factorisation; dividing by the greatest common divisor reduces any fraction to it.
MathWorld — Sylvester's Four-Point Problem. The probability that four points chosen independently and uniformly in a region have a quadrilateral as convex hull is 25/36 for the square.
Wikipedia — Bentley–Ottmann algorithm. J. Bentley and T. Ottmann (1979) report the k crossings of n segments in O((n + k) log n) time, an improvement on testing all pairs when k = o(n²/log n).
C++ — Euler165.cpp. The linked immutable revision contains cross, intersection_point, count_true_intersections and run_checkpoints; it is the implementation record for this article.

Problem 165 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <set>
#include <string>
#include <tuple>
#include <vector>

namespace {

using i64 = long long;

struct Options {
    int segments = 5000;
    bool run_checkpoints = true;
};

struct Point {
    i64 x = 0;
    i64 y = 0;
};

struct Segment {
    Point a;
    Point b;
};

i64 cross(const Point& a, const Point& b, const Point& c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

int sgn(const i64 x) {
    return (x > 0) - (x < 0);
}

struct RationalPoint {
    i64 x_num = 0;
    i64 x_den = 1;
    i64 y_num = 0;
    i64 y_den = 1;

    bool operator<(const RationalPoint& other) const {
        return std::tie(x_num, x_den, y_num, y_den) < std::tie(other.x_num, other.x_den, other.y_num, other.y_den);
    }
};

RationalPoint intersection_point(const Segment& s1, const Segment& s2) {
    const i64 a1 = s1.b.y - s1.a.y;
    const i64 b1 = s1.a.x - s1.b.x;
    const i64 c1 = a1 * s1.a.x + b1 * s1.a.y;

    const i64 a2 = s2.b.y - s2.a.y;
    const i64 b2 = s2.a.x - s2.b.x;
    const i64 c2 = a2 * s2.a.x + b2 * s2.a.y;

    i64 den = a1 * b2 - a2 * b1;
    i64 x_num = c1 * b2 - c2 * b1;
    i64 y_num = a1 * c2 - a2 * c1;

    if (den < 0) {
        den = -den;
        x_num = -x_num;
        y_num = -y_num;
    }

    const i64 gx = std::gcd(std::llabs(x_num), den);
    const i64 gy = std::gcd(std::llabs(y_num), den);

    RationalPoint out;
    out.x_num = x_num / gx;
    out.x_den = den / gx;
    out.y_num = y_num / gy;
    out.y_den = den / gy;
    return out;
}

std::vector<Segment> generate_segments(const int n) {
    std::vector<i64> s(static_cast<std::size_t>(4 * n + 1), 0);
    s[0] = 290797;
    for (int i = 0; i < 4 * n; ++i) {
        s[static_cast<std::size_t>(i + 1)] = (s[static_cast<std::size_t>(i)] * s[static_cast<std::size_t>(i)]) % 50515093;
    }

    std::vector<Segment> segs(static_cast<std::size_t>(n));
    for (int i = 0; i < n; ++i) {
        segs[static_cast<std::size_t>(i)].a.x = s[static_cast<std::size_t>(4 * i + 1)] % 500;
        segs[static_cast<std::size_t>(i)].a.y = s[static_cast<std::size_t>(4 * i + 2)] % 500;
        segs[static_cast<std::size_t>(i)].b.x = s[static_cast<std::size_t>(4 * i + 3)] % 500;
        segs[static_cast<std::size_t>(i)].b.y = s[static_cast<std::size_t>(4 * i + 4)] % 500;
    }
    return segs;
}

i64 count_true_intersections(const std::vector<Segment>& segs) {
    std::set<RationalPoint> points;
    const int segments = static_cast<int>(segs.size());

    for (int i = 0; i < segments; ++i) {
        for (int j = i + 1; j < segments; ++j) {
            const i64 c1 = cross(segs[static_cast<std::size_t>(i)].a, segs[static_cast<std::size_t>(i)].b, segs[static_cast<std::size_t>(j)].a);
            const i64 c2 = cross(segs[static_cast<std::size_t>(i)].a, segs[static_cast<std::size_t>(i)].b, segs[static_cast<std::size_t>(j)].b);
            const i64 c3 = cross(segs[static_cast<std::size_t>(j)].a, segs[static_cast<std::size_t>(j)].b, segs[static_cast<std::size_t>(i)].a);
            const i64 c4 = cross(segs[static_cast<std::size_t>(j)].a, segs[static_cast<std::size_t>(j)].b, segs[static_cast<std::size_t>(i)].b);

            if (sgn(c1) * sgn(c2) >= 0 || sgn(c3) * sgn(c4) >= 0) {
                continue;
            }

            points.insert(intersection_point(segs[static_cast<std::size_t>(i)], segs[static_cast<std::size_t>(j)]));
        }
    }

    return static_cast<i64>(points.size());
}

i64 solve(const int segments) {
    const std::vector<Segment> segs = generate_segments(segments);
    return count_true_intersections(segs);
}

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--segments=", options.segments)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.segments >= 2;
}

bool run_checkpoints() {
    const std::vector<Segment> sample = generate_segments(1);
    if (sample[0].a.x != 27 || sample[0].a.y != 144 || sample[0].b.x != 12 || sample[0].b.y != 232) {
        std::cerr << "Checkpoint failed for BBS seed" << '\n';
        return false;
    }

    std::vector<Segment> stated(3);
    stated[0] = Segment{Point{27, 44}, Point{12, 32}};
    stated[1] = Segment{Point{46, 53}, Point{17, 62}};
    stated[2] = Segment{Point{46, 70}, Point{22, 40}};
    if (count_true_intersections(stated) != 1LL) {
        std::cerr << "Checkpoint failed for stated sample segments" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.segments) << '\n';
    return 0;
}

Python

from math import gcd

def solve():
    num_segments = 5000
    s = [0] * (4 * num_segments + 1)
    s[0] = 290797
    for i in range(4 * num_segments):
        s[i + 1] = (s[i] * s[i]) % 50515093

    segs = []
    for i in range(num_segments):
        ax = s[4 * i + 1] % 500
        ay = s[4 * i + 2] % 500
        bx = s[4 * i + 3] % 500
        by = s[4 * i + 4] % 500
        segs.append((ax, ay, bx, by))

    def cross(ax, ay, bx, by, cx, cy):
        return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)

    def sgn(x):
        return (x > 0) - (x < 0)

    def intersection_point(s1, s2):
        a1 = s1[3] - s1[1]
        b1 = s1[0] - s1[2]
        c1 = a1 * s1[0] + b1 * s1[1]
        a2 = s2[3] - s2[1]
        b2 = s2[0] - s2[2]
        c2 = a2 * s2[0] + b2 * s2[1]
        den = a1 * b2 - a2 * b1
        x_num = c1 * b2 - c2 * b1
        y_num = a1 * c2 - a2 * c1
        if den < 0:
            den = -den
            x_num = -x_num
            y_num = -y_num
        gx = gcd(abs(x_num), den)
        gy = gcd(abs(y_num), den)
        return (x_num // gx, den // gx, y_num // gy, den // gy)

    points = set()
    for i in range(len(segs)):
        si = segs[i]
        for j in range(i + 1, len(segs)):
            sj = segs[j]
            c1 = cross(si[0], si[1], si[2], si[3], sj[0], sj[1])
            c2 = cross(si[0], si[1], si[2], si[3], sj[2], sj[3])
            c3 = cross(sj[0], sj[1], sj[2], sj[3], si[0], si[1])
            c4 = cross(sj[0], sj[1], sj[2], sj[3], si[2], si[3])
            if sgn(c1) * sgn(c2) >= 0 or sgn(c3) * sgn(c4) >= 0:
                continue
            points.add(intersection_point(si, sj))

    return str(len(points))

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler165 {
    public static void main(String[] args) {
        int[] s = new int[20001];
        s[0] = 290797;
        for (int i = 1; i <= 20000; i++)
            s[i] = (int) ((long) s[i - 1] * s[i - 1] % 50515093);
        int[] t = new int[20000];
        for (int i = 0; i < 20000; i++)
            t[i] = s[i + 1] % 500;
        int[][] segs = new int[5000][4];
        for (int n = 0; n < 5000; n++) {
            segs[n][0] = t[4 * n];
            segs[n][1] = t[4 * n + 1];
            segs[n][2] = t[4 * n + 2];
            segs[n][3] = t[4 * n + 3];
        }
        Set<Long> pts = new HashSet<>();
        for (int i = 0; i < 5000; i++) {
            long x1 = segs[i][0], y1 = segs[i][1], x2 = segs[i][2], y2 = segs[i][3];
            long dx1 = x2 - x1, dy1 = y2 - y1;
            for (int j = i + 1; j < 5000; j++) {
                long x3 = segs[j][0], y3 = segs[j][1], x4 = segs[j][2], y4 = segs[j][3];
                long dx2 = x4 - x3, dy2 = y4 - y3;
                long denom = dx1 * dy2 - dy1 * dx2;
                if (denom == 0)
                    continue;
                long tn = (x3 - x1) * dy2 - (y3 - y1) * dx2;
                long un = (x3 - x1) * dy1 - (y3 - y1) * dx1;
                if (denom > 0) {
                    if (tn <= 0 || tn >= denom || un <= 0 || un >= denom)
                        continue;
                } else {
                    if (tn >= 0 || tn <= denom || un >= 0 || un <= denom)
                        continue;
                }
                // Exact point: reduce fraction
                long pxN = x1 * denom + tn * dx1, pxD = denom;
                long pyN = y1 * denom + tn * dy1, pyD = denom;
                long g1 = gcd(Math.abs(pxN), Math.abs(pxD));
                pxN /= g1;
                pxD /= g1;
                long g2 = gcd(Math.abs(pyN), Math.abs(pyD));
                pyN /= g2;
                pyD /= g2;
                if (pxD < 0) {
                    pxN = -pxN;
                    pxD = -pxD;
                }
                if (pyD < 0) {
                    pyN = -pyN;
                    pyD = -pyD;
                }
                // Encode as unique key
                pts.add(hash4(pxN, pxD, pyN, pyD));
            }
        }
        System.out.println(pts.size());
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static long hash4(long a, long b, long c, long d) {
        return ((a * 1000003 + b) * 1000003 + c) * 1000003 + d;
    }
}