-- Constructing the unique circumcenter of a noncollinear triangle.
-- Part of Cumulus, the Grace17 mathematics library.
-- Copyright (c) 2026 Grant Passmore.
-- Released under the Apache License 2.0. See LICENSE in the repository root.
-- Authors: Grant Passmore

module Grace.Cumulus.Geometry.Circumcenter
import TriangleCenters
import LinearAlgebra.TwoByTwo

namespace Euclidean
export open Euclidean

-- Right-hand side of the linear perpendicular-bisector equation.
def bisectorOffset (a b : Point) : Real :=
  half(x(b) * x(b) + y(b) * y(b) - (x(a) * x(a) + y(a) * y(a)))

theorem bisectorOffset.real (a b : Point)
  shows Real(bisectorOffset(a, b))
proof
  unfold bisectorOffset
  typecheck

theorem bisectorOffset.double (a b : Point)
  shows bisectorOffset(a, b) + bisectorOffset(a, b) =
    x(b) * x(b) + y(b) * y(b) - (x(a) * x(a) + y(a) * y(a))
proof
  have doubled := half.add_self at x(b) * x(b) + y(b) * y(b) - (x(a) * x(a) + y(a) * y(a)) with *
  simp [bisectorOffset]

-- Squared-distance equality is exactly one linear equation in the center.
theorem circumcenter.bisector_equation (a b o : Point)
  shows Euclidean.distanceSquared(o, a) = Euclidean.distanceSquared(o, b) ↔
    (x(b) - x(a)) * x(o) + (y(b) - y(a)) * y(o) = bisectorOffset(a, b)
proof
  have offsetReal := bisectorOffset.real at a, b with *
  have doubled := bisectorOffset.double at a, b with *
  have firstRadius := Euclidean.distanceSquared.coordinates at o, a with *
  have secondRadius := Euclidean.distanceSquared.coordinates at o, b with *
  constructor
  . intro radii
    rw [firstRadius, secondRadius] at radii
    apply realAddValue.double_cancel with offsetReal
    linear_combination radii - doubled
  . intro equation
    rw [firstRadius, secondRadius]
    linear_combination equation + equation + doubled

-- Cramer's rule is total; correctness and uniqueness require noncollinearity.
def circumcenter (a b c : Set) : Set :=
  (LinearSystem2.solutionX(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)),
   LinearSystem2.solutionY(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)))

theorem circumcenter.point (a b c : Point)
  shows Point(circumcenter(a, b, c))
proof
  unfold circumcenter
  exact pair.point with *

theorem circumcenter.coordinates (a b c : Point)
  shows x(circumcenter(a, b, c)) = LinearSystem2.solutionX(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)) ∧
    y(circumcenter(a, b, c)) = LinearSystem2.solutionY(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c))
proof
  have horizontalProjection := x.pair at LinearSystem2.solutionX(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)), LinearSystem2.solutionY(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)) with *
  have verticalProjection := y.pair at LinearSystem2.solutionX(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)), LinearSystem2.solutionY(x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c)) with *
  constructor
  . simp [circumcenter]
  . simp [circumcenter]

theorem circumcenter.correct (a b c : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  shows Circumcenter(circumcenter(a, b, c), a, b, c)
proof
  have nonsingular : (x(b) - x(a)) * (y(c) - y(a)) - (y(b) - y(a)) * (x(c) - x(a)) ≠ 0
    intro zero
    apply noncollinear
    simp [Collinear, directionCross]
  have equations := LinearSystem2.correct at x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c) with *
  have point := circumcenter.point at a, b, c with *
  have coordinates := circumcenter.coordinates at a, b, c with *
  rcases coordinates with ⟨horizontal, vertical⟩
  rcases equations with ⟨first, second⟩
  rw [← horizontal, ← vertical] at first
  rw [← horizontal, ← vertical] at second
  have firstRadii := circumcenter.bisector_equation.backward at a, b, circumcenter(a, b, c) with *
  have secondRadii := circumcenter.bisector_equation.backward at a, c, circumcenter(a, b, c) with *
  unfold Circumcenter
  refine ⟨firstRadii, secondRadii⟩

theorem circumcenter.unique (a b c o p : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  requires firstCenter : Circumcenter(o, a, b, c)
  requires secondCenter : Circumcenter(p, a, b, c)
  shows o = p
proof
  have nonsingular : (x(b) - x(a)) * (y(c) - y(a)) - (y(b) - y(a)) * (x(c) - x(a)) ≠ 0
    intro zero
    apply noncollinear
    simp [Collinear, directionCross]
  rcases firstCenter with ⟨firstRadii, secondRadii⟩
  rcases secondCenter with ⟨otherFirstRadii, otherSecondRadii⟩
  have first := circumcenter.bisector_equation.forward at a, b, o with *
  have second := circumcenter.bisector_equation.forward at a, c, o with *
  have otherFirst := circumcenter.bisector_equation.forward at a, b, p with *
  have otherSecond := circumcenter.bisector_equation.forward at a, c, p with *
  have equalCoordinates := LinearSystem2.unique at x(b) - x(a), y(b) - y(a), x(c) - x(a), y(c) - y(a), bisectorOffset(a, b), bisectorOffset(a, c), x(o), y(o), x(p), y(p) with *
  rcases equalCoordinates with ⟨horizontal, vertical⟩
  exact point.ext with *

theorem circumcenter.exists (a b c : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  shows ∃ o, Point(o) ∧ Circumcenter(o, a, b, c)
proof
  have point := circumcenter.point at a, b, c with *
  have circle := circumcenter.correct with *
  refine ⟨circumcenter(a, b, c), point, circle⟩

-- No center is supplied by the caller: all three are constructed foundationally.
theorem euler_line.exists (a b c : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  shows ∃ o, ∃ g, ∃ h, Point(o) ∧ Point(g) ∧ Point(h) ∧
    Circumcenter(o, a, b, c) ∧ Centroid(g, a, b, c) ∧ Orthocenter(h, a, b, c) ∧
    Collinear(h, g, o) ∧ MedianRatio(h, g, o)
proof
  have oPoint := circumcenter.point at a, b, c with *
  have circle := circumcenter.correct with *
  have gPoint := centroid.point at a, b, c with *
  have center := centroid.correct at a, b, c with *
  have line := euler_line at circumcenter(a, b, c), a, b, c, centroid(a, b, c) with *
  rcases line with ⟨h, hPoint, altitudes, collinear, ratio⟩
  refine ⟨circumcenter(a, b, c), centroid(a, b, c), h, oPoint, gPoint, hPoint, circle, center, altitudes, collinear, ratio⟩

-- Every noncollinear triangle has a center, constructed from its circumcenter.
theorem orthocenter.exists (a b c : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  shows ∃ h, Point(h) ∧ Orthocenter(h, a, b, c)
proof
  have oPoint := circumcenter.point at a, b, c with *
  have circle := circumcenter.correct with *
  have hPoint := orthocenter.point at circumcenter(a, b, c), a, b, c with *
  have altitudes := orthocenter.correct with *
  refine ⟨orthocenter(circumcenter(a, b, c), a, b, c), hPoint, altitudes⟩

-- An intersection of two altitudes lies on the third as well.
theorem orthocenter.of_two_altitudes (a b c h : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  requires firstAltitude : PerpendicularDirections(a, h, b, c)
  requires secondAltitude : PerpendicularDirections(b, h, c, a)
  shows Orthocenter(h, a, b, c)
proof
  have center := orthocenter.exists with *
  rcases center with ⟨k, kPoint, altitudes⟩
  obtain ⟨otherFirstAltitude, otherSecondAltitude, otherThirdAltitude⟩ := altitudes
  have same := altitudes.unique at a, b, c, h, k with *
  rwa same

-- The explicit construction agrees with every geometric orthocenter.
theorem orthocenter.eq_constructed (a b c h : Point)
  requires noncollinear : ¬Collinear(a, b, c)
  requires center : Orthocenter(h, a, b, c)
  shows h = orthocenter(circumcenter(a, b, c), a, b, c)
proof
  have oPoint := circumcenter.point at a, b, c with *
  have circle := circumcenter.correct with *
  have hPoint := orthocenter.point at circumcenter(a, b, c), a, b, c with *
  have altitudes := orthocenter.correct with *
  exact orthocenter.unique with *

end Euclidean