Projective geometry

Stand between the rails of a straight railway line in the Karoo and look along it. The rails, which never meet, appear to converge to a point on the horizon. The sleepers, laid exactly one metre apart, appear closer and closer together. A broadcast camera above a football pitch in Buenos Aires does the same thing to the touchlines and the penalty box: lengths, angles and parallelism are all distorted.

Yet some things survive every photograph. Straight lines stay straight. Points on a line stay on a line, and lines through a point still pass through a point. And a less obvious quantity survives too: from four points on a line one can compute a number, the cross-ratio, that is the same on the ground and in the picture. Sports-analytics software uses exactly these invariants to recover the true positions of players from television images, and photogrammetry uses them to measure buildings from photographs.

Projective geometry is the geometry of what perspective preserves. To make it work, we must give parallel lines a place to meet: a point at infinity for each direction, all lying on a line at infinity — the horizon. The reward is a geometry of remarkable symmetry, in which every theorem comes with a free twin, and in which ellipses, parabolas and hyperbolas turn out to be the same curve seen from different places.

Put your eye at the origin of $\R^3$ and hold a transparent picture plane at height $z = 1$. Each point $(x, y)$ of the picture is seen along a line through the eye, the line through $(x, y, 1)$. So points of the picture are lines through the origin. A line in the picture, together with the eye, spans a plane through the origin, so lines of the picture are planes through the origin.

Now the lines through the eye that are horizontal — parallel to the picture plane — never meet it. They are perfectly good lines through the origin, so we add them as new points: the points at infinity, one for each horizontal direction, and together they form a new line, the horizontal plane through the eye, which is the line at infinity. Two parallel lines in the picture span two planes through the eye whose intersection is a horizontal line: they now meet, at a point at infinity, exactly as the rails meet on the horizon.

This is the real projective plane $\R P^2$: the set of lines through the origin in $\R^3$. Linear algebra on $\R^3$ becomes geometry on $\R P^2$. Two distinct lines through the origin span a plane: two points determine a line. Two distinct planes through the origin meet in a line: two lines determine a point, with no exception for parallels. The two statements are twins, and that symmetry runs through the whole subject.

Throughout, $K$ is a field, usually $\R$.

The last sentence of the proof is the principle of duality. Incidence, $\ell \cdot p = 0$, is symmetric in $\ell$ and $p$. So exchanging the words "point" and "line" (and "collinear" and "concurrent", "join" and "meet") turns every true statement about incidence in $KP^2$ into another true statement. Theorem 1 is self-dual; Desargues' theorem below is dual to its own converse, and Pascal's theorem on conics is dual to Brianchon's. We get every projective theorem twice for the price of one proof.

On the projective line the same argument shows that a projective transformation is determined by the images of three distinct points. On $\R P^2$ the four-point theorem says that any convex quadrilateral can be mapped onto a square by a perspective map — the homography used to rectify photographs.

Three points and a cross-ratio determine the fourth point, which is how one locates the fourth fence post in a photograph from the first three (first experiment). Permuting the four points changes the cross-ratio in a controlled way: the $24$ orderings give at most six values,

which collapse to three when $\lambda = -1$. Four points with cross-ratio $-1$ form a harmonic range: for example $0$ and $\infty$ separate $1$ and $-1$ harmonically, and the midpoint $m$ of $a$ and $b$ satisfies $(a, b; m, \infty) = -1$. This is how projective geometry recovers midpoints once the line at infinity is known — exactly what a perspective artist does with the horizon.

There is a classical synthetic proof too. If the two triangles lie in different planes of a three-dimensional projective space, the lines $AB$ and $A'B'$ lie in the plane $OAB$ and therefore meet, and all three meeting points lie on the line where the planes $ABC$ and $A'B'C'$ intersect. The planar case follows by projecting a three-dimensional configuration. So Desargues' theorem in the plane is a shadow of the third dimension.

The proof used that multiplication commutes (in forming $\alpha\alpha'$ and writing the cross product). That is no accident. Projective planes can be built over any division ring $D$, where multiplication need not commute, such as Hamilton's quaternions. Hilbert (1899) showed that Desargues' theorem holds in a projective plane exactly when it can be coordinatised by a division ring, and Pappus' theorem holds exactly when that division ring is commutative, a field. So Desargues holds and Pappus fails in the quaternionic plane. There are even planes satisfying the incidence axioms in which Desargues fails, such as Moulton's plane (1902) and finite planes of order $9$. Hessenberg (1905) proved synthetically that Pappus implies Desargues, and by Wedderburn's theorem (1905) every finite division ring is a field, so in a finite plane Desargues implies Pappus.

For example $x^2 + y^2 - z^2$ becomes $y^2 - uw$ under $x = \frac{u - w}2$, $z = \frac{u + w}2$: in the chart $w = 1$ that is the parabola $y^2 = u$. So the circle and the parabola are the same projective curve, seen with different lines chosen as "at infinity". The visualisation below moves the line at infinity continuously. Two further classical facts, stated without proof: five points in general position determine exactly one conic (a symmetric $3 \times 3$ matrix up to scale has five degrees of freedom, and each point imposes one linear condition); and Pascal's theorem (1640) says that for a hexagon inscribed in a conic, the three meets of opposite sides are collinear. Pappus' theorem is the case of the degenerate conic $\ell \cup m$.

Mathlib's Projectivization K V (notation ℙ K V) is the type of lines through $0$ in $V$, with Projectivization.mk K v hv the point with homogeneous coordinates $v \ne 0$. Homogeneous coordinates are defined up to scale, and opposite vectors give the same point — one point at infinity per direction:

/-- **Homogeneous coordinates are defined up to scale**: `[v] = [w]` in `ℙ K V` exactly when
`v = a • w` for some scalar `a` (Mathlib's `Projectivization.mk_eq_mk_iff'`). -/
theorem mk_eq_mk_iff {K V : Type*} [Field K] [AddCommGroup V] [Module K V] (v w : V)
    (hv : v ≠ 0) (hw : w ≠ 0) :
    Projectivization.mk K v hv = Projectivization.mk K w hw ↔ ∃ a : K, a • w = v :=
  Projectivization.mk_eq_mk_iff' K v w hv hw

/-- **Opposite directions give the same point.** `[v] = [−v]`: in particular a direction in the
plane and its reverse determine the *same* point at infinity `[x : y : 0] = [−x : −y : 0]`. -/
theorem mk_neg {K : Type*} [Field K] (v : Fin 3 → K) (hv : v ≠ 0) :
    Projectivization.mk K v hv = Projectivization.mk K (-v) (neg_ne_zero.mpr hv) := by
  rw [Projectivization.mk_eq_mk_iff]
  exact ⟨-1, by simp⟩

Theorem 1 with the cross product $u ⨯_3 v$ from Mathlib.LinearAlgebra.CrossProduct. The theorem line_unique is proved once, and read in two ways: for points $p, q$ it says the line is unique; for lines it says the meet is unique. That is the principle of duality inside a proof assistant.

/-- **The line through `p` and `q` is `p × q`**: the cross product is orthogonal to both factors,
so both points lie on it (Mathlib's `dot_self_cross` and `dot_cross_self`). Read dually, with
`p`, `q` lines: the point `p × q` lies on both lines. -/
theorem cross_on_both {K : Type*} [CommRing K] (p q : Fin 3 → K) :
    p ⬝ᵥ (p ⨯₃ q) = 0 ∧ q ⬝ᵥ (p ⨯₃ q) = 0 :=
  ⟨dot_self_cross p q, dot_cross_self p q⟩

/-- `p × q ≠ 0` exactly when `p` and `q` are linearly independent, i.e. represent **distinct**
points (Mathlib's `crossProduct_ne_zero_iff_linearIndependent`). -/
theorem cross_ne_zero_iff {K : Type*} [Field K] (p q : Fin 3 → K) :
    p ⨯₃ q ≠ 0 ↔ LinearIndependent K ![p, q] :=
  crossProduct_ne_zero_iff_linearIndependent

/-- If the line `ℓ` passes through `p` and `q`, then `ℓ × (p × q) = 0`, i.e. `ℓ` and `p × q` are
parallel vectors. Proof: the expansion `ℓ × (p × q) = (ℓ · q) p − (p · ℓ) q` (Mathlib's
`cross_cross_eq_smul_sub_smul'`) has both coefficients zero. -/
theorem incident_cross_zero {K : Type*} [CommRing K] (ℓ p q : Fin 3 → K) (hp : ℓ ⬝ᵥ p = 0)
    (hq : ℓ ⬝ᵥ q = 0) : ℓ ⨯₃ (p ⨯₃ q) = 0 := by
  rw [cross_cross_eq_smul_sub_smul', hq, dotProduct_comm, hp, zero_smul, zero_smul, sub_zero]

/-- **Two distinct points lie on exactly one line.** If `p × q ≠ 0` (distinct points) and `ℓ`
passes through both, then `ℓ` is a scalar multiple of `p × q` — the same projective line. By
duality (swap the roles of points and lines) the *same theorem* says that two distinct lines
meet in exactly one point, their cross product. -/
theorem line_unique {K : Type*} [Field K] (ℓ p q : Fin 3 → K) (hpq : p ⨯₃ q ≠ 0)
    (hp : ℓ ⬝ᵥ p = 0) (hq : ℓ ⬝ᵥ q = 0) : ∃ c : K, c • (p ⨯₃ q) = ℓ := by
  have h1 : (p ⨯₃ q) ⨯₃ ℓ = 0 := by
    rw [← cross_anticomm, incident_cross_zero ℓ p q hp hq, neg_zero]
  have h2 : ¬ LinearIndependent K ![p ⨯₃ q, ℓ] := by
    rw [← crossProduct_ne_zero_iff_linearIndependent, not_not]
    exact h1
  rw [LinearIndependent.pair_iff' hpq] at h2
  push Not at h2
  exact h2

/-- **Collinearity test**: three points `u, v, w` of `ℙ²` are collinear (lie on a common line)
exactly when `det [u, v, w] = 0`; the triple product `u · (v × w)` is that determinant
(Mathlib's `triple_product_eq_det`). -/
theorem triple_product_det {K : Type*} [CommRing K] (u v w : Fin 3 → K) :
    u ⬝ᵥ (v ⨯₃ w) = Matrix.det ![u, v, w] :=
  triple_product_eq_det u v w

The cross-ratio and Theorem 3, in the affine form of the proof (a Möbius map and the cancellation of $(ad - bc)^2$):

/-- The **cross-ratio** of four numbers (affine coordinates of four points on a line):
`(a, b; c, d) = ((a − c)(b − d)) / ((a − d)(b − c))`. (Lean's division by `0` gives `0`.) -/
def crossRatio {K : Type*} [Field K] (a b c d : K) : K :=
  ((a - c) * (b - d)) / ((a - d) * (b - c))

/-- The key identity behind cross-ratio invariance: for `f(x) = (a x + b)/(c x + d)`,
`f(x) − f(y) = (ad − bc)(x − y) / ((c x + d)(c y + d))`. -/
theorem mobius_sub {K : Type*} [Field K] (a b c d x y : K) (hx : c * x + d ≠ 0)
    (hy : c * y + d ≠ 0) :
    (a * x + b) / (c * x + d) - (a * y + b) / (c * y + d)
      = (a * d - b * c) * (x - y) / ((c * x + d) * (c * y + d)) := by
  rw [div_sub_div _ _ hx hy]
  congr 1
  ring

/-- **The cross-ratio is a projective invariant.** A projective transformation of the line,
`f(x) = (a x + b)/(c x + d)` with `ad − bc ≠ 0`, preserves the cross-ratio of any four points
where it is defined. In each of the four differences the factor `ad − bc` and the denominators
`c x + d` cancel between numerator and denominator — exactly the lesson's proof. -/
theorem cross_ratio_mobius {K : Type*} [Field K] (a b c d p q r s : K) (hdet : a * d - b * c ≠ 0)
    (hp : c * p + d ≠ 0) (hq : c * q + d ≠ 0) (hr : c * r + d ≠ 0) (hs : c * s + d ≠ 0) :
    let f := fun x => (a * x + b) / (c * x + d)
    crossRatio (f p) (f q) (f r) (f s) = crossRatio p q r s := by
  intro f
  simp only [crossRatio, f]
  rw [mobius_sub a b c d p r hp hr, mobius_sub a b c d q s hq hs, mobius_sub a b c d p s hp hs,
    mobius_sub a b c d q r hq hr]
  set D := a * d - b * c
  have hP : (c * p + d) * (c * r + d) * ((c * q + d) * (c * s + d)) ≠ 0 := by
    apply mul_ne_zero <;> apply mul_ne_zero <;> assumption
  rw [div_mul_div_comm, div_mul_div_comm]
  rw [show (c * p + d) * (c * s + d) * ((c * q + d) * (c * r + d))
      = (c * p + d) * (c * r + d) * ((c * q + d) * (c * s + d)) by ring]
  rw [div_div_div_cancel_right₀ hP]
  rw [show D * (p - r) * (D * (q - s)) = (D * D) * ((p - r) * (q - s)) by ring,
    show D * (p - s) * (D * (q - r)) = (D * D) * ((p - s) * (q - r)) by ring]
  exact mul_div_mul_left _ _ (mul_ne_zero hdet hdet)

Desargues' theorem exactly as proved above, and Pappus' theorem in homogeneous coordinates adapted to the two lines. The Pappus statement covers points at infinity too (it allows any $[a_1 : 0 : a_2]$ on the first line), and its proof is a single polynomial identity checked by ring. Both hold in any commutative ring, and changing one point off its line makes the ring check fail.

/-- **Desargues' theorem** in homogeneous coordinates. Triangles `ABC` and `A'B'C'` are in
perspective from `O`: `A' = O + a A`, `B' = O + b B`, `C' = O + c C` (every point of the line `OA`
other than `A` has a representative of this form). Then `P = aA − bB` lies on the lines `AB` and
`A'B'` (because `P = A' − B'`), similarly `Q` on `BC` and `B'C'` and `R` on `CA` and `C'A'`, and
`P + Q + R = 0` makes `P, Q, R` collinear: the triple product `P · (Q × R)` vanishes. -/
theorem desargues {K : Type*} [CommRing K] (O A B C : Fin 3 → K) (a b c : K) :
    let A' := O + a • A
    let B' := O + b • B
    let C' := O + c • C
    let P := a • A - b • B
    let Q := b • B - c • C
    let R := c • C - a • A
    (A ⬝ᵥ (B ⨯₃ P) = 0 ∧ A' ⬝ᵥ (B' ⨯₃ P) = 0) ∧
    (B ⬝ᵥ (C ⨯₃ Q) = 0 ∧ B' ⬝ᵥ (C' ⨯₃ Q) = 0) ∧
    (C ⬝ᵥ (A ⨯₃ R) = 0 ∧ C' ⬝ᵥ (A' ⨯₃ R) = 0) ∧
    P ⬝ᵥ (Q ⨯₃ R) = 0 := by
  intro A' B' C' P Q R
  simp only [A', B', C', P, Q, R, cross_apply, dotProduct, Fin.sum_univ_three]
  refine ⟨⟨?_, ?_⟩, ⟨?_, ?_⟩, ⟨?_, ?_⟩, ?_⟩ <;> simp <;> ring

/-- **Pappus' theorem** in homogeneous coordinates. After a projective change of coordinates the
two lines are `y = 0` and `x = 0`, so `A, B, C = [a₁ : 0 : a₂], …` and `A', B', C' = [0 : d₁ : d₂], …`.
The three cross-joins `AB' ∩ A'B`, `AC' ∩ A'C`, `BC' ∩ B'C` — each a cross product of two
cross products — are collinear: their triple product is identically zero. The proof is a
polynomial identity checked by `ring`, valid in every *commutative* ring: commutativity of
multiplication is exactly what Pappus needs. -/
theorem pappus {K : Type*} [CommRing K] (a₁ a₂ b₁ b₂ c₁ c₂ d₁ d₂ e₁ e₂ f₁ f₂ : K) :
    let A : Fin 3 → K := ![a₁, 0, a₂]
    let B : Fin 3 → K := ![b₁, 0, b₂]
    let C : Fin 3 → K := ![c₁, 0, c₂]
    let A' : Fin 3 → K := ![0, d₁, d₂]
    let B' : Fin 3 → K := ![0, e₁, e₂]
    let C' : Fin 3 → K := ![0, f₁, f₂]
    let P := (A ⨯₃ B') ⨯₃ (A' ⨯₃ B)
    let Q := (A ⨯₃ C') ⨯₃ (A' ⨯₃ C)
    let R := (B ⨯₃ C') ⨯₃ (B' ⨯₃ C)
    P ⬝ᵥ (Q ⨯₃ R) = 0 := by
  intro A B C A' B' C' P Q R
  simp only [A, B, C, A', B', C', P, Q, R, cross_apply, dotProduct, Fin.sum_univ_three]
  simp
  ring

Conics: the explicit change of coordinates taking the circle to a parabola, and Mathlib's version of Sylvester's law used in Theorem 6(a):

/-- **The circle and the parabola are projectively the same conic.** With
`P = [[1/2, 0, −1/2], [0, 1, 0], [1/2, 0, 1/2]]` (invertible: `det P = 1/2`), the form
`x² + y² − z²` (the circle `x² + y² = 1` when `z = 1`) becomes `y² − u w`, i.e. `Pᵀ C P = D` —
the parabola `y² = u` in the chart `w = 1`. Concretely, `x = (u − w)/2`, `z = (u + w)/2`. -/
theorem circle_to_parabola :
    let P : Matrix (Fin 3) (Fin 3) ℚ := !![1/2, 0, -1/2; 0, 1, 0; 1/2, 0, 1/2]
    let C : Matrix (Fin 3) (Fin 3) ℚ := !![1, 0, 0; 0, 1, 0; 0, 0, -1]
    let D : Matrix (Fin 3) (Fin 3) ℚ := !![0, 0, -1/2; 0, 1, 0; -1/2, 0, 0]
    Pᵀ * C * P = D ∧ P.det ≠ 0 := by
  intro P C D
  refine ⟨?_, ?_⟩
  · ext i j
    fin_cases i <;> fin_cases j <;> norm_num [P, C, D, Matrix.mul_apply, Fin.sum_univ_three]
  · norm_num [P, Matrix.det_fin_three]

/-- **Sylvester's classification of real quadratic forms**: a non-degenerate quadratic form on a
real `n`-dimensional space is equivalent to `Σ wᵢ xᵢ²` with every `wᵢ = ±1` (Mathlib's
`QuadraticForm.equivalent_one_neg_one_weighted_sum_squared`). For `n = 3`, up to an overall sign
only `x² + y² + z²` (an empty conic) and `x² + y² − z²` (the circle) remain: so all non-empty
non-degenerate real conics are projectively equivalent. -/
theorem real_quadratic_forms {M : Type*} [AddCommGroup M] [Module ℝ M] [FiniteDimensional ℝ M]
    (Q : QuadraticForm ℝ M) (hQ : (QuadraticMap.associated (R := ℝ) Q).SeparatingLeft) :
    ∃ w : Fin (Module.finrank ℝ M) → ℝ,
      (∀ i, w i = -1 ∨ w i = 1) ∧ QuadraticMap.Equivalent Q (QuadraticMap.weightedSumSquares ℝ w) :=
  Q.equivalent_one_neg_one_weighted_sum_squared hQ

What is and is not formalised. Independent: $[v] = [-v]$, the uniqueness half of Theorem 1 (and its dual), the cross-ratio identity and Theorem 3 for Möbius maps, Desargues' theorem (Theorem 4) and Pappus' theorem (Theorem 5) as coordinate identities, and the circle-to-parabola transformation. From Mathlib: equality in ℙ K V, the cross-product incidence facts, the triple product as a determinant, and the classification of real quadratic forms. Not formalised: Theorem 2 (four points), the reduction "choose coordinates so that the lines are the axes" (we state Pappus directly in those coordinates), Theorem 6(b), Hilbert's coordinatisation theorems and Pascal's theorem. Mathlib does not yet have synthetic projective planes. The Lean exercises ask for $[v] = [-v]$ and the incidence lemma.

The first program stores points of the projective line as homogeneous pairs of fractions, so that $\infty = [1 : 0]$ needs no special case. It applies $300$ random Möbius transformations and checks that every cross-ratio is preserved, while ordinary ratios of lengths almost never are. It lists the six values of the cross-ratio under reordering (three in the harmonic case), and locates the fourth fence post of a photograph from the first three.

"""The cross-ratio: invariant under every projective map of the line.

Points of the projective line are stored as homogeneous pairs [x : y] of exact
fractions, so the point at infinity [1 : 0] needs no special treatment. With
the bracket [p q] = p_x q_y - p_y q_x, the cross-ratio is

    (a, b; c, d) = [a c][b d] / ([a d][b c]),

which agrees with ((a - c)(b - d)) / ((a - d)(b - c)) for finite points. The
program applies seeded random Moebius transformations x -> (ax + b)/(cx + d)
(2 x 2 matrices with non-zero determinant) to random quadruples and checks
the cross-ratio exactly, while the simple ratio (a - c)/(b - c) changes. It
lists the six values taken by the cross-ratio under the 24 orderings, and uses
invariance to locate a fence post in a perspective photograph.

All arithmetic is exact: a finite exhaustive check of the cases listed.
"""

import random
from fractions import Fraction
from itertools import permutations

INFINITY = (Fraction(1), Fraction(0))

def point(x):
    """The homogeneous pair [x : 1] of a finite point."""
    return (Fraction(x), Fraction(1))

def bracket(p, q):
    """The 2 x 2 determinant [p q] = p_x q_y - p_y q_x."""
    return p[0] * q[1] - p[1] * q[0]

def cross_ratio(a, b, c, d):
    """(a, b; c, d) = [a c][b d] / ([a d][b c]) for homogeneous pairs."""
    return (bracket(a, c) * bracket(b, d)) / (bracket(a, d) * bracket(b, c))

def apply(m, p):
    """Apply the 2 x 2 matrix m = ((a, b), (c, d)) to the homogeneous pair p."""
    (a, b), (c, d) = m
    return (a * p[0] + b * p[1], c * p[0] + d * p[1])

def affine(p):
    """The affine coordinate x/y of [x : y], or None for the point at infinity."""
    return None if p[1] == 0 else p[0] / p[1]

def random_mobius(rng):
    """A random integer 2 x 2 matrix with non-zero determinant."""
    while True:
        m = tuple(tuple(Fraction(rng.randint(-9, 9)) for _ in range(2))
                  for _ in range(2))
        if m[0][0] * m[1][1] - m[0][1] * m[1][0] != 0:
            return m

def six_values(lam):
    """The values of the cross-ratio over all 24 orderings of four points.

    The points lam, 1, 0, infinity have cross-ratio (lam, 1; 0, inf) = lam.
    """
    base = [point(lam), point(1), point(0), INFINITY]
    values = {cross_ratio(*[base[i] for i in perm])
              for perm in permutations(range(4))}
    return sorted(values)

def fourth_point(a, b, c, lam):
    """The affine x with ((a - c)(b - x)) / ((a - x)(b - c)) = lam.

    Clearing denominators gives b - x = k (a - x) with k = lam (b - c)/(a - c).
    """
    k = lam * (b - c) / (a - c)
    return (k * a - b) / (k - 1)

def photo(x):
    """Position in a perspective photo of a point x metres along the fence."""
    return Fraction(x) / (x + 2)

def main():
    rng = random.Random(2024)
    trials, invariant, ratio_changed = 300, 0, 0
    for _ in range(trials):
        xs = rng.sample(range(-20, 21), 4)
        quad = [point(x) for x in xs]
        if rng.random() < 0.25:
            quad[rng.randrange(4)] = INFINITY
        m = random_mobius(rng)
        images = [apply(m, p) for p in quad]
        if any(bracket(images[i], images[j]) == 0
               for i in range(4) for j in range(i)):
            continue
        invariant += cross_ratio(*images) == cross_ratio(*quad)
        a, b, c = (affine(p) for p in quad[:3])
        a2, b2, c2 = (affine(p) for p in images[:3])
        if None not in (a, b, c, a2, b2, c2):
            ratio_changed += (a - c) / (b - c) != (a2 - c2) / (b2 - c2)
    print(f"{trials} random Moebius maps on random quadruples (exact):")
    print(f"   cross-ratio unchanged in {invariant} of {trials} cases")
    print(f"   simple ratio (a - c)/(b - c) changed in {ratio_changed} cases")
    print()
    for lam in (Fraction(3), Fraction(-1)):
        vals = six_values(lam)
        print(f"orderings of four points with cross-ratio {lam}: "
              f"{len(vals)} values")
        print("   " + ", ".join(str(v) for v in vals))
    print()
    posts = [0, 1, 2, 3]
    cr_ground = cross_ratio(*[point(x) for x in posts])
    seen = [photo(x) for x in posts[:3]]
    print("Fence posts 1 m apart at 0, 1, 2, 3 have cross-ratio", cr_ground)
    print("In the photo (x -> x/(x + 2)) the first three appear at",
          ", ".join(str(s) for s in seen))
    x = fourth_point(*seen, cr_ground)
    print(f"Invariance predicts the fourth at {x}; the photo has"
          f" {photo(posts[3])}")
    print()
    print("Epistemic status: exact rational arithmetic over finitely many\n"
          "random cases (and two worked examples). This is evidence, not a\n"
          "proof, of invariance; the proof is the algebra in the lesson,\n"
          "checked in Lean as cross_ratio_mobius.")

if __name__ == "__main__":
    main()

The second program builds the finite planes $\mathbb F_p P^2$ for $p = 2, 3, 5, 7$. It checks by exhaustion that each has $p^2 + p + 1$ points and lines, that two points determine one line and two lines one point, and that the conic $x^2 + y^2 = z^2$ has exactly $p + 1$ points for odd $p$, just as a real conic is a circle, the projective line in disguise. It then tests Desargues and Pappus on random configurations in $\mathbb F_7P^2$.

"""Finite projective planes PG(2, p): counting, incidence, Desargues, Pappus.

Over the field F_p of integers mod a prime p, the projective plane PG(2, p)
has as points the non-zero vectors of F_p^3 up to non-zero scalars, and lines
likewise (a line l contains the point x when l . x = 0 mod p). The program

* counts points and lines (p^2 + p + 1 each) and points per line (p + 1);
* checks exhaustively that two distinct points lie on exactly one line and
  two distinct lines meet in exactly one point;
* counts the points of the conic x^2 + y^2 = z^2 (p + 1 for odd p);
* tests the theorems of Desargues and Pappus on seeded random configurations
  in PG(2, 7), computing joins and meets with cross products mod p.

Incidence checks are exhaustive for the listed p; the theorem checks are
random samples, i.e. evidence for that plane only.
"""

import random
from itertools import combinations, product

def normalise(v, p):
    """Scale a non-zero vector mod p so that its first non-zero entry is 1."""
    v = [x % p for x in v]
    for x in v:
        if x:
            inv = pow(x, p - 2, p)
            return tuple(y * inv % p for y in v)
    raise ValueError("zero vector")

def points(p):
    """All points of PG(2, p), as normalised triples."""
    return sorted({normalise(v, p) for v in product(range(p), repeat=3)
                   if any(v)})

def dot(u, v, p):
    return sum(a * b for a, b in zip(u, v)) % p

def cross(u, v, p):
    """The cross product mod p: the line through two points (or the meet)."""
    return ((u[1] * v[2] - u[2] * v[1]) % p, (u[2] * v[0] - u[0] * v[2]) % p,
            (u[0] * v[1] - u[1] * v[0]) % p)

def collinear(u, v, w, p):
    """True when the three points lie on one line (determinant zero mod p)."""
    return dot(u, cross(v, w, p), p) == 0

def incidence_report(p):
    """(#points, set of points-per-line counts, axioms hold) for PG(2, p)."""
    pts = points(p)
    lines = pts  # by duality the lines are the same triples
    on = {l: [x for x in pts if dot(l, x, p) == 0] for l in lines}
    per_line = {len(v) for v in on.values()}
    two_points = all(sum(1 for l in lines if dot(l, x, p) == 0 and
                         dot(l, y, p) == 0) == 1
                     for x, y in combinations(pts, 2))
    two_lines = all(sum(1 for x in pts if dot(l, x, p) == 0 and
                        dot(m, x, p) == 0) == 1
                    for l, m in combinations(lines, 2))
    return len(pts), per_line, two_points and two_lines

def conic_points(p):
    """Number of points of x^2 + y^2 = z^2 in PG(2, p)."""
    return sum(1 for (x, y, z) in points(p) if (x * x + y * y - z * z) % p == 0)

def desargues_trial(rng, p):
    """Random triangles perspective from a point; are the meets collinear?"""
    while True:
        o, a, b, c = (rng.choices(range(p), k=3) for _ in range(4))
        s, t, u = (rng.randrange(1, p) for _ in range(3))
        a2 = [(o[i] + s * a[i]) % p for i in range(3)]
        b2 = [(o[i] + t * b[i]) % p for i in range(3)]
        c2 = [(o[i] + u * c[i]) % p for i in range(3)]
        tri = [a, b, c, a2, b2, c2, o]
        if any(not any(v) for v in tri):
            continue
        if collinear(a, b, c, p) or collinear(a2, b2, c2, p):
            continue
        meets = [cross(cross(a, b, p), cross(a2, b2, p), p),
                 cross(cross(b, c, p), cross(b2, c2, p), p),
                 cross(cross(c, a, p), cross(c2, a2, p), p)]
        if any(not any(m) for m in meets):
            continue
        return collinear(*meets, p)

def pappus_trial(rng, p):
    """Three points on each of two lines; are the cross-joins collinear?"""
    while True:
        l, m = (rng.choices(range(p), k=3) for _ in range(2))
        if not any(cross(l, m, p)):
            continue
        pts = points(p)
        on_l = [x for x in pts if dot(l, x, p) == 0 and dot(m, x, p) != 0]
        on_m = [x for x in pts if dot(m, x, p) == 0 and dot(l, x, p) != 0]
        a, b, c = rng.sample(on_l, 3)
        a2, b2, c2 = rng.sample(on_m, 3)
        meets = [cross(cross(a, b2, p), cross(a2, b, p), p),
                 cross(cross(a, c2, p), cross(a2, c, p), p),
                 cross(cross(b, c2, p), cross(b2, c, p), p)]
        return collinear(*meets, p)

def main():
    print(" p   points  points/line  2 pts -> 1 line, 2 lines -> 1 pt"
          "  conic pts")
    for p in (2, 3, 5, 7):
        n, per_line, axioms = incidence_report(p)
        conic = conic_points(p) if p > 2 else "-"
        print(f"{p:2d}   {n:6d}  {sorted(per_line)!s:>11}  {axioms!s:>31}"
              f"  {conic!s:>9}")
    print("(p^2 + p + 1 points, p + 1 per line; a conic has p + 1 points)")
    print()
    rng = random.Random(2024)
    trials = 300
    des = sum(desargues_trial(rng, 7) for _ in range(trials))
    pap = sum(pappus_trial(rng, 7) for _ in range(trials))
    print(f"PG(2, 7): Desargues held in {des} of {trials} random "
          "configurations")
    print(f"PG(2, 7): Pappus held in {pap} of {trials} random configurations")
    print()
    print("Epistemic status: the counts and incidence axioms are exhaustive\n"
          "checks for p = 2, 3, 5, 7. The Desargues and Pappus results are\n"
          "random samples in one finite plane: evidence, not proof. The\n"
          "proofs (valid over every field) are in the lesson and in Lean.")

if __name__ == "__main__":
    main()

The incidence counts are exhaustive, hence proofs for those particular planes. The random Desargues and Pappus tests are evidence in one plane; the coordinate proofs (and the Lean theorems) cover every field at once.

The curve below is the image of the unit circle under the projective transformation $[x : y : z] \mapsto [x : y : z + ay]$, drawn in the chart $z = 1$. That is, the circle is fixed and the line sent to infinity moves: it is the line $y = -1/a$. For $0 \le a < 1$ that line misses the circle and the image is an ellipse. At $a = 1$ the line is tangent at $(0, -1)$, and the image is a parabola. For $a > 1$ it cuts the circle twice, and the image is a hyperbola whose two branches are joined through the two points at infinity. Slide $a$ slowly through $1$.

Möbius transformations $z \mapsto \frac{az + b}{cz + d}$ are the projective transformations of the complex projective line $\C P^1 = \C \cup \{\infty\}$, the Riemann sphere. They preserve the complex cross-ratio, and they send circles and lines to circles and lines (a line is a circle through $\infty$). Watch the grid lines bend into circles.