Circle–Line Intersection — Mathematical Foundations and Numerical Analysis¶

Objects.
  Circle  Γ = { P ∈ ℝ² : |P − C|² = r² },  C ∈ ℝ², r > 0.
  Line    ℓ = { O + t·D : t ∈ ℝ },          O ∈ ℝ², D ∈ ℝ², D ≠ 0.

Intersection set.
  Γ ∩ ℓ = { O + t·D : |O + t·D − C|² = r² }.

Let f = O − C. Substituting P = O + t·D:
  |t·D + f|² = r²
  (D·D) t² + 2 (f·D) t + (f·f − r²) = 0.                         (★)

Quadratic form a·t² + b·t + c = 0 with
  a = D·D  > 0,
  b = 2 (f·D),
  c = f·f − r².

Discriminant  Δ = b² − 4ac.

Cardinality.
  |Γ ∩ ℓ| = 0  ⇔ Δ < 0     (line misses Γ),
  |Γ ∩ ℓ| = 1  ⇔ Δ = 0     (ℓ tangent to Γ),
  |Γ ∩ ℓ| = 2  ⇔ Δ > 0     (ℓ secant to Γ).

Geometric quantities (perpendicular method).
  d  = dist(C, ℓ) = |f − ((f·D)/(D·D)) D|   (length of the component of f ⟂ D),
  F  = O + ((f·D)/(D·D)) D... wait, foot uses (C−O): F = O + ((C−O)·D / D·D) D,
  h  = √(r² − d²)   (real only when d ≤ r).

Claim. For a non-degenerate line (D ≠ 0) and circle (r > 0), the algebraic
roots t = (−b ± √Δ)/(2a) of (★) yield exactly the points F ± h·D̂ of the
geometric method, and Δ = 4a²(r² − d²).

Setup. WLOG translate so C = 0 (translation preserves both methods). Then
f = O, and (★) reads (D·D)t² + 2(O·D)t + (O·O − r²) = 0.

Decompose O along and perpendicular to D. Let D̂ = D/|D|, |D| = √a. Write
  O = α D̂ + β n̂,
where n̂ ⟂ D̂, |n̂| = 1. Then:
  O·D  = α |D| = α √a,         (since D̂·D = |D|, n̂·D = 0)
  O·O  = α² + β².
The perpendicular distance from C = 0 to ℓ is exactly d = |β|, the component
of O orthogonal to the line. (β is the signed offset.)

Compute Δ.
  b  = 2(O·D)   = 2α√a
  c  = O·O − r² = α² + β² − r²
  Δ  = b² − 4ac = 4aα² − 4a(α² + β² − r²)
                = 4a(r² − β²)
                = 4a(r² − d²).                                   (†)

Hence sign(Δ) = sign(r² − d²): the count tests agree. ✔

Also (†) gives, for a normalised direction (a = 1):  Δ = 4(r² − d²).

Roots. The midpoint of the two roots is t* = −b/(2a) = −α√a/a = −α/√a.
The point at t* is
  P(t*) = O + t* D = α D̂ + β n̂ + (−α/√a)(√a D̂) = α D̂ + β n̂ − α D̂ = β n̂.
That is the foot of perpendicular F (the closest point on ℓ to C = 0):
its component along D̂ is 0, so it is where the perpendicular from C meets ℓ. ✔

Half-width in t. The roots are t* ± √Δ/(2a). Convert the t-offset to a length
along the line (multiply by |D| = √a):
  Δt · |D| = (√Δ/(2a))·√a = √Δ/(2√a) = √(4a(r²−d²)) / (2√a)
           = (2√a · √(r²−d²)) / (2√a) = √(r² − d²) = h.            (‡)

So each algebraic root sits a distance h on either side of F along the line —
precisely F ± h·D̂, the geometric construction. The two methods are identical. ∎

Δ = b² − 4ac,  with a = D·D > 0.

Sensitivity. How does the count flip as inputs vary?
  ∂Δ/∂(r²) = −4a·∂c/∂(r²)... since c = f·f − r², ∂c/∂(r²) = −1, so ∂Δ/∂(r²) = 4a > 0.
  → growing the circle (↑ r²) increases Δ monotonically: a miss becomes a
    tangent then a secant. The tangent radius is r_tan = d (Δ = 0 ⇔ r = d).

  ∂Δ/∂d² = −4a < 0 (from Δ = 4a(r² − d²) with a = D·D):
  → pushing the line away from the centre (↑ d) decreases Δ: secant → tangent
    → miss. The tangent occurs exactly at d = r.

Degenerate boundaries.
  a → 0  (D → 0): (★) ceases to be quadratic; the "line" collapses to a point.
                  Handle separately (degeneracy taxonomy).
  r → 0:  c = f·f ≥ 0, Δ = b² − 4a(f·f) = −4a·d²·... ≤ 0; intersection only if
          d = 0 (point C lies on the line) → single point. Consistent with a
          zero-radius circle being the point C.

Conditioning of the count. The decision |Γ∩ℓ| ∈ {0,1,2} is governed by the
sign of Δ. Near Δ = 0 the sign is ill-conditioned: a perturbation δ in inputs
shifts Δ by O(‖input‖·δ), and when |Δ| is comparable to that shift the sign is
unreliable. This is *exactly* the geometric statement that "is this line tangent?"
is numerically hard — there is no way around it; one can only choose a principled
tolerance band (see error bounds).

The textbook formula t = (−b ± √Δ)/(2a) is numerically poor. When b and √Δ
have the same magnitude and the same sign, the root using "−b + √Δ" (or
"−b − √Δ") computes a small number as the difference of two near-equal large
numbers → catastrophic cancellation, losing most significant digits.

Stable formulation (Numerical Recipes). Compute one root with addition (no
cancellation), the other via Vieta's product t1·t2 = c/a.

  q = −½ ( b + sign(b)·√Δ )          # sign(b) makes |q| as large as possible
  t1 = q / a
  t2 = c / q                          # since t1·t2 = c/a ⇒ t2 = (c/a)/t1 = c/q

Both roots are now computed without subtracting near-equal quantities.
If b = 0 the formula degenerates gracefully to t = ±√(−c/a) (= ±√(r²−d²)/√a · …),
which is already well-conditioned.

Halved-b refinement. Using b' = f·D (so b = 2b'), Δ' = b'² − a·c and
  q = −( b' + sign(b')·√Δ' );  t1 = q/a;  t2 = c/q.
Fewer multiplications, identical stability.

Model. IEEE-754 double; unit round-off u ≈ 1.11e−16. Each ⊕,⊖,⊗ introduces a
relative error ≤ u. We bound the error in Δ' = b'² − a·c.

Computed Δ̂' = fl(fl(b'·b') − fl(a·c)). Standard analysis gives
  |Δ̂' − Δ'| ≤ u·(|b'|² + |a c|) + O(u²)
           ≲ u·(|b'|² + |a||c|).

When the true Δ' is small (near tangent) but |b'|² and |a c| are large and
nearly equal, the *relative* error in Δ' blows up:
  rel_err(Δ') ≈ u · (|b'|² + |a c|) / |Δ'|  →  ∞  as Δ' → 0.

This is unavoidable: the condition number of "subtract two near-equal numbers"
is (|x|+|y|)/|x−y|. No reformulation removes it; it is intrinsic to deciding
tangency. Consequences and the only real mitigations:

1. SCALE REDUCTION. Translate C to the origin (f = O − C with small |O|) and
   normalise D (a = 1). Then |b'|² + |a c| is minimised for the given geometry,
   shrinking the absolute error floor. This is the highest-leverage fix.

2. RELATIVE TOLERANCE. Declare tangency when
      |Δ'| ≤ τ · (|b'|² + |a c|),     τ a few × u  (e.g. 1e−10 .. 1e−12).
   This makes the band proportional to the achievable precision, not absolute.

3. HIGHER PRECISION where it matters. Compute Δ' in extended precision
   (e.g. Kahan/two-product + two-sum, or float128) for the count decision only;
   the points themselves tolerate double precision once the count is fixed.

4. POINT ERROR. Given a correctly-classified secant, the points P = O + t·D have
   error ≈ u·(|O| + |t|·|D|) per coordinate — well-conditioned away from tangency.
   The fragility is entirely in the *count*, not in the *coordinates*.

Concrete demonstration of why the naive root formula loses digits, and how the
c/q form rescues it. Take a grazing ray nearly tangent to a unit-normalised line.

Setup (already in local frame, a = D·D = 1):
  a  = 1
  b  = 2.000000000          (so b' = 1.0)
  c  = 0.99999999999999     (chosen so Δ is tiny)
  Δ  = b² − 4ac = 4 − 3.99999999999996 = 4.0e−14   (true value)
  √Δ ≈ 2.0e−7

Naive roots t = (−b ± √Δ)/(2a):
  t_small = (−2 + 2.0e−7)/2 = (−1.9999998)/2 = −0.99999990    (large, fine)
  t_large = (−2 − 2.0e−7)/2 = −1.00000010                     (large, fine)
  ... here BOTH roots are ≈ −1, near each other; the DIFFERENCE between them,
  which is what physically matters (the chord half-width), is 2.0e−7 computed as
  the difference of two numbers each ≈ 1.0 — only ~7 of 16 digits survive.

Stable form q = −(b + sign(b)·√Δ)/2 = −(2 + 2.0e−7)/2 = −1.00000010:
  t1 = q/a = −1.00000010
  t2 = c/q = 0.99999999999999 / (−1.00000010) = −0.99999990
  The product t1·t2 = c/a is preserved EXACTLY by construction, so the small
  separation is recovered to full precision instead of being differenced away.

Lesson: the failure is not in computing each root's MAGNITUDE (both ≈ 1) but in
recovering their tiny SEPARATION. Vieta's product (t2 = c/q) sidesteps the
subtraction entirely. This is the whole reason production solvers never use the
textbook ± formula directly for both roots.

The perpendicular distance d ties this topic to the orientation/cross-product
primitive of 02-line-intersection, unifying the geometry toolkit.

For the line through A, B and point C:
  cross(A,B,C) = (Bx−Ax)(Cy−Ay) − (By−Ay)(Cx−Ax)  = 2·SignedArea(△ABC).

The triangle's area is also ½·base·height = ½·|AB|·d. Equating:
  |cross(A,B,C)| = |AB|·d    ⇒    d = |cross(A,B,C)| / |AB|.

So the SAME scalar that gives "left/right of the line" (orientation) gives the
DISTANCE to the line once divided by the base length. The sign of cross(A,B,C)
additionally tells you WHICH SIDE of the line the centre lies on — useful for
oriented queries (e.g. is the circle entering from the left or right half-plane).

Consequence for exact predicates: because cross(A,B,C) is an exact integer for
integer inputs, the quantity |AB|²·d² = cross² is exact, and comparing it to
|AB|²·r² decides the count with NO division and NO sqrt:

  count test (exact):   sign( |AB|²·r²  −  cross(A,B,C)² )
                      = sign( a·r²·... )  ≡  sign(Δ')   up to the positive factor a.

This is the cleanest exact formulation: one integer cross product, squared,
compared against r² scaled by the squared base length.

When inputs are integers (or rationals), the COUNT can be decided EXACTLY,
avoiding floating-point entirely. With integer C, r, and integer line points
A, B:

  Let D = B − A, f = A − C   (integer vectors).
  a = D·D, b' = f·D, c = f·f − r²   (all exact integers; use wide types).
  Δ' = b'² − a·c   (exact integer).

  count = 0 if Δ' < 0;  1 if Δ' = 0;  2 if Δ' > 0.   ← EXACT, no rounding.

This is the robust-predicate analogue of the orientation test from
02-line-intersection: the *decision* is an integer sign test. Overflow is the
only enemy — products of coordinates near 10^k reach 10^(4k) in Δ', so use
__int128 / BigInteger / Python big ints when coordinates are large.

The POINTS, however, are generally irrational (they involve √Δ'), so reporting
exact coordinates requires an algebraic-number representation (a + b√Δ'). In
practice: decide the count exactly with integers, then compute the point
coordinates in double precision — the count was the only ill-conditioned part.

Adaptive precision (Shewchuk-style). Compute Δ' in double first with a forward
error bound; if |Δ̂'| exceeds the bound, the sign is certified — return it.
Only when |Δ̂'| is within the error bound do you escalate to exact integer /
extended-precision evaluation. This pays the exact-arithmetic cost only on the
rare near-tangent inputs, keeping the common case fast.

A second derivation, useful when the line arrives as ax + by + c0 = 0 (implicit
/ normal form) rather than two points. Normalise so a² + b² = 1; then (a, b) is
the unit normal and |a·x + b·y + c0| is exactly the perpendicular distance from
(x, y) to the line.

Distance from centre C = (cx, cy):
  d = |a·cx + b·cy + c0|          (since a² + b² = 1).

Foot of perpendicular:
  F = C − (a·cx + b·cy + c0)·(a, b).     (step from C along the unit normal)

Intersection points (when d ≤ r):
  Let û = (−b, a)   (unit direction ALONG the line, ⟂ to the normal).
  h = √(r² − d²).
  P = F ± h·û.

Equivalence to the parametric derivation is immediate: the parametric direction
D, normalised, equals ±(−b, a), and the foot/half-chord construction is identical.
The implicit form simply hands you the unit normal directly, so no projection
dot-product is needed — one fewer step when the line is already in normal form.

This is the form most convenient inside half-plane / clipping pipelines
(12-half-plane-intersection), where lines are naturally stored as (a, b, c0).

Beyond the COUNT (sign of Δ), how accurate are the POINT COORDINATES? We bound
the first-order sensitivity of a root t to perturbations.

From a·t² + b·t + c = 0, implicit differentiation w.r.t. a parameter p gives
  (2a·t + b)·(∂t/∂p) + t²·(∂a/∂p) + t·(∂b/∂p) + (∂c/∂p) = 0.
Note 2a·t + b = ±√Δ (the derivative of the quadratic at a root). Hence
  ∂t/∂p = −( t²·∂a/∂p + t·∂b/∂p + ∂c/∂p ) / (±√Δ).

The denominator √Δ → 0 at tangency: the roots become INFINITELY SENSITIVE to
input perturbation exactly at the tangent point. This is the analytic statement
of the same fact the error analysis gave numerically — near tangency, both the
count AND the precise location of the (merging) points are ill-conditioned.

Away from tangency (Δ bounded away from 0), ∂t/∂p is O(1/√Δ) and well-behaved;
the point coordinates inherit the input precision. So:
  • secant, well clear of tangent:  points accurate to ~unit round-off.
  • approaching tangent:             location accuracy degrades like 1/√Δ.
  • at tangent:                      a single point, but its position is the
                                     most sensitive case — small input noise
                                     moves it most.