Hart’s theorem

Hart’s theorem says

If a triangle be formed by the arcs of three circles, the inscribed and the three escribed circles are all tangent to a new circle or line.

Here “triangle” means a three-sided figure whose sides are portions of a circle. The inscribed circle is the largest circle that can fit inside the three-sided figure.

The “escribed” circles are analogous to the excircles in the previous post: you extend two sides and find a circle that is tangent to the triangle side and the extended side. The difference here being that the side extensions are now circles.

Incircles and Excircles of Pythagorean triangles

This post will reveal the connection between my two previous posts: one on the Star Trek lemma and one on Pythagorean triples.

In the process of writing the latter, I looked at the Wikipedia article on Pythagorean triples and noticed this curious paragraph.

In every Pythagorean triangle, the radius of the incircle and the radii of the three excircles are positive integers. Specifically, for a primitive triple the radius of the incircle is r = n(m − n), and the radii of the excircles opposite the sides m2 − n22mn, and the hypotenuse m2 + n2 are respectively m(m − n), n(m + n), and m(m + n).

The citation for the paragraph above was the book by my former officemate, which led to the post on the Star Trek lemma. The passage in Arthur Baragar’s book that Wikipedia cites is Exercise 15.3.

Let ΔABC be a right angle triangle with sides of integer length. Prove that the inradius r and the exradii ra, rb, and rc are all integers.

I don’t know whether Arthur discovered this theorem, but I’ll call it Baragar’s theorem for this post.

Incircles and excircles

To unpack Baragar’s theorem, let’s start by saying what incircles and excircles are. Incircles are more familiar. The incircle of a triangle is the largest circle that can be inscribed inside the triangle, and the radius of this circle is the inradius.

Since an incircle is an inscribed circle, you might expect an excircle to be a circumscribed circle, but that’s not it. There are three excircles, one for each side. To find the excircle for a side, extend the other two sides and find the circle tangent to the side and the two extensions. The radius of an excircle is its exradius.

Proof

Baragar’s theorem follows directly from Euclid’s formula for Pythagorean triples mentioned in the previous post

\begin{align*} a &= m^2 - n^2 \\ b &= 2mn \\ c &= m^2 + n^2 \end{align*}

and formulas for the inradius r and the exradii ra, rb, and rc.

\begin{align*} r &= \frac{K}{s} \\ r_a &= \frac{K}{s - a} \\ r_b &= \frac{K}{s - b} \\ r_c &= \frac{K}{s - c} \\ \end{align*}

Here K is the area of the triangle, which in our case is ab/2, and s is the semiperimeter, half the perimeter.

Expressing the radii in terms of m and n gives the values cited by Wikipedia above.

Illustrating the theorem

I’d like to write a Python script to illustrate the theorem, and knowing the radii of the circles help, but we also need to know the centers of the circles.

The center of the incircle is the weighted average of the vertices, with weights given by the lengths of the opposite sides. That is, if the vertices are AB, and C, and the sides opposite these vertices are ab, and c, the the incenter is

I = \frac{aA + bB + cC}{a + b + c}

The centers for the excircles have remarkably similar expressions. For the incenter of the circle opposite a vertex, flip the sign of the corresponding side.

\begin{align*} I_a = \frac{-aA + bB + cC}{-a + b + c} \\ I_a = \frac{aA - bB + cC}{a - b + c} \\ I_a = \frac{aA + bB - cC}{a + b - c} \\ \end{align*}

Python code

Putting it all together, here’s an illustrate the theorem.

And here’s the code that produced it. Note that everything in this section works for right triangles in general, not just Pythagorean triangles.

import numpy as np
import matplotlib.pyplot as plt

def connect(A, B):
    plt.plot([A[0], B[0]], [A[1], B[1]], "C0")

def draw_circle(c, r, color):
    t = np.linspace(0, 2*np.pi)
    plt.plot(r*np.cos(t) + c[0], r*np.sin(t) + c[1], color=color)

a, b, c, = 3, 4, 5

A = np.array([0, b])
B = np.array([-a, 0])
C = np.array([0, 0])

s = (a + b + c)/2
K = a*b/2
r = K/s
ra = K/(s - a)
rb = K/(s - b)
rc = K/(s - c)
I = (a*A + b*B + c*C)/(a + b + c)
Ia = (-a*A + b*B + c*C)/(-a + b + c)
Ib = (a*A - b*B + c*C)/(a - b + c)
Ic = (a*A + b*B - c*C)/(a + b - c)

draw_circle(I, r, "C1")
draw_circle(Ia, ra, "C2")
draw_circle(Ib, rb, "C3")
draw_circle(Ic, rc, "C4")

plt.plot([-2*rc, 2*rb], [0, 0], "C0")
plt.plot([0, 0], [-2*ra, 2*rc], "C0")
plt.plot([(-2*ra - b)*a/b, 2*rb], [-2*ra, 2*rb*b/a + b], "C0")

plt.gca().set_aspect("equal")
plt.show()

Consecutive Pythagorean triangle sides

In this post we find all Pythagorean triples that contain consecutive numbers, all Pythagorean triples (abc) such that a + 1 = b or b + 1 = c.

a + 1 = b

George Osborne wrote a paper [1] addressing the question of when the squares of two consecutive numbers is also a square. Geometrically this is asking for primitive Pythagorean triples for which the legs are consecutive integers.

He proved that the sequence shorter legs satisfies the recurrence relation

u_{n+2} = 6 u_{n+1} - u_{n+1} + 2
with initial conditions u0 = 0 and u1 = 1. This is OEIS sequence A001652.

The method for solving recurrences like the one above is analogous to the method for solving linear differential equations. See a solution here. This gives us the following formula for the terms:

u_n = \dfrac{1 + \sqrt{2}}{4} \left(3 + 2\sqrt{2}\right)^n + \dfrac{1 - \sqrt{2}}{4} \left(3 - 2\sqrt{2}\right)^n - \dfrac{1}{2}

b + 1 = c

It’s also possible for the longer side and hypotenuse of a Pythagorean triangle to be consecutive numbers, as in the (5, 12, 13) triangle.

All primitive Pythagorean triples are given by Euclid’s formula

\begin{align*} a &= m^2 - n^2 \\ b &= 2mn \\ c &= m^2 + n^2 \end{align*}

with integers m > n > 0. If b and c are consecutive numbers, then

c - b = 1 = m^2 + n^2 - 2mn = (m -n)^2

and so mn + 1. Therefore all possible values of b are given by 2n(n + 1) for n > 1.

[1] Geo. A. Osborne. A Problem in Number Theory. The American Mathematical Monthly, Vol. 21, No. 5 (May, 1914), pp. 148-150

The Star Trek lemma

I was reading an article this evening and saw a footnote to a book by Arthur Baragar [1]. This caught my eye because he was my officemate at UT for a year.

I found his book on Archive.org and was surprised to see “The Star Trek Lemma” in the table of contents. What could this be?

It’s a theorem that goes back to Euclid that applies to an angle formed by connecting a point to two other points on a circle. The theorem says “The measure of an inscribed angle is half the measure of the arc it subtends.” But why call it the Star Trek lemma? Quoting Arthur:

In the spirit of Euclid, we will refer to this theorem as the Star Trek lemma because of the figure associated with the statement of the theorem. … Before Star Trek, as far as I know, this theorem had no name, though some might call it Euclid III.20, which is its proposition number in Euclid’s Elements (Book III, Proposition 20).

Here is my reconstruction of the figure given in the book.

Baragar's illustration of the Star Trek lemma

The lemma says that ∠BAC is half of ∠BOC.

[1] Baragar, Arthur (2001), A Survey of Classical and Modern Geometries: With Computer Activities, Prentice Hall

 

Regular expressions that work “everywhere”

The most frustrating aspect of regular expressions is that implementations vary. Features supported in one tool may not be supported at all in another tool, or they may be supported with slightly different syntax.

I learned regular expressions in the context Perl, a maximalist regex environment. This led to frustration when features I expect to work are missing [1]. One way around this is to use Perl analogs of other tools, but this is very non-standard. I want to be able to send colleagues and clients code that works out of the box.

As I mentioned in my post on computational survivalism, I occasionally need to work on computers that I cannot install software on. So a better approach is to identify a subset of regex features that work everywhere. The stricter your definition of “everywhere” the less this includes. The strictest subset would be

  • literals
  • character classes […]
  • the special characters . * ^ $

A more relaxed definition of “everywhere” would be the tools you most care about. Currently the tools I most want to use with regular expressions are sed, awk, grep, and Emacs.

Awk as lowest common denominator

If you use the Gnu versions of sed, awk, and grep, and use the -E option with sed and grep, then the list of common features is bigger. The regular expression features of of the three tools are similar, and awk’s features are supported in the other tools, with one exception: word boundaries in awk are \< and \> rather than \b and \B.

I wrote about Awk’s regex features here.

Emacs as the oddball

Emacs supports analogs of most of awk’s regex features. However, the characters

    + ? ( ) { } |

all require a backslash in front in order to act like the awk counterparts. Also, the analog of \s and \S in awk is \s- and \S- in Emacs.

Instead of meaning space or nonspace, \s and \S in Emacs begin a (negated) character class, and one of those classes is - for space. But there are many others. For example, \s. stands for a punctuation character and \S. stands for a non-punctuation character.

What works everywhere

So for my definition of “everywhere,” with the caveats mentioned above, the following features work everywhere. YMMV.

    .
    ^, $
    […], [^…]
    *
    \w, \W, \s, \S
    \1 - \9 backreferences
    \b \B
    ? + 
    | alternation
    {n,m} for counting matches
    (...) capturing

One footnote is that gawk supports backreferences in replacement strings but not in regular expressions per se.

[1] To some extent, basic Perl features work elsewhere and advanced features do not, depending on your idea of what is basic or advanced. I think of look-around features as advanced, and that tracks. But I think of \d for digits as basic, but that’s not supported in many regex flavors.

Lobachevsky’s integral formula

Let f be an even function with period π. Then the following remarkable theorem by Lobachevsky holds.

\int_0^\infty \frac{\sin^2 x}{x^2} f(x) \, dx = \int_0^\infty\frac{\sin x} x f(x) \, dx = \int_0^{\pi/2} f(x) \, dx

This theorem is useful in Fourier analysis and signal processing. It’s useful to know even in the special case f(x) = 1.

For a “jinc” analog, see this paper.

***

Every time I see the name Lobachevsky I think of Tom Lehrer’s song about him. You can find the words here and the audio here.

Related posts

Queens on a prime order board

The n queens problem is to place on an n × n chessboard n queens so that none attacks any other. This means there is only one queen on every horizontal, vertical, and diagonal line.

When n is a prime number ≥ 5, it is sufficient to place the queens on a line that has slope 2, 3, 4, …, n − 2. (The slope cannot be 1 because that’s a diagonal. And it cannot be n − 1 because n − 1 = −1 mod n is also a diagonal.) [1]

Here we imagine opposite edges of the board being joined together. Geometrically, this makes the chessboard a torus (donut). Algebraically, the points on a line of slope s have the coordinates

(akbks)

where addition is carried out mod n.

All solutions to the n queens problem have this form when n = 5. Some solutions will have this form for larger prime values of n but not all.

For example, when n = 7, here is a solution where all the queens are on a line of slope 2.

But here is another solution where the queens do not all lie on a line of constant slope.

Related posts

[1] W. H. Bussey. A Note on the Problem of the Eight Queens. The American Mathematical Monthly, Vol. 29, No. 7 (August 1922), pp. 252–253

All pieces on a 6 by 5 board

I’ve written a couple posts lately on getting an LLM to generate code to solve chess problems. The first used Claude to generate Prolog and the second used ChatGPT to generate Prolog. This post will use Claude to generate Z3/Python code.

The puzzle is one I’ve written about before:

Place all the pieces—king, queen, two bishops, two knights, and two rooks—on a 6 × 5 chessboard, with the requirement that the two bishops be on opposite colored squares and no piece is attacking another.

Incidentally, it’s common for “piece” to exclude pawns, as above. But then what do you call all the things on a chessboard? You might call them “chess pieces,” in which case a pawn is a “chess piece” but not a “piece.” One convention is to use “chessmen” or simply “men” to include pieces and pawns.

This was the prompt I used.

Write Z3/Python code to find all solutions to the following chess puzzle.

Place all the pieces—king, queen, two bishops, two knights, and two rooks—on a 6 × 5 chessboard, with the requirement that the two bishops be on opposite colored squares and no piece is attacking another.

The code found 192 = 8 × 24 solutions. The factor of 8 comes from 23 ways of swapping the pairs of bishops, knights, and rooks. The script reports

Total raw solutions: 192
Unique solutions (deduplicating piece-pair swaps): 24

── Solution 1 ──
  0 1 2 3 4
0 . K . . N
1 . . . . B
2 . . R . .
3 Q . . . .
4 . . . R .
5 . B . . N
  King: (0,1)
  Queen: (3,0)
  Bishop1: (5,1) [light]
  Bishop2: (1,4) [dark]
  Knight1: (5,4)
  Knight2: (0,4)
  Rook1: (4,3)
  Rook2: (2,2)

── Solution 2 ──
  0 1 2 3 4
0 N . . K .
1 B . . . .
2 . . R . .
3 . . . . Q
4 . R . . .
5 N . . B .
  King: (0,3)
  Queen: (3,4)
  Bishop1: (5,3) [light]
  Bishop2: (1,0) [dark]
  Knight1: (5,0)
  Knight2: (0,0)
  Rook1: (4,1)
  Rook2: (2,2)

...

── Solution 24 ──
  0 1 2 3 4
0 B . . B .
1 N . . . .
2 . . . . Q
3 . . R . .
4 . R . . .
5 N . . K .
  King: (5,3)
  Queen: (2,4)
  Bishop1: (0,0) [light]
  Bishop2: (0,3) [dark]
  Knight1: (5,0)
  Knight2: (1,0)
  Rook1: (3,2)
  Rook2: (4,1)

Here’s the code.

"""
Chess Puzzle Solver using Z3
Place king, queen, 2 bishops, 2 knights, 2 rooks on a 6×5 board
such that:
  - The two bishops are on opposite-colored squares
  - No piece attacks any other piece
"""

from z3 import (
    Int, Solver, Or, And, Not, Distinct, If,
    sat, ForAll, Implies
)
import itertools

ROWS = 6
COLS = 5
TOTAL_CELLS = ROWS * COLS  # 30

PIECES = ["king", "queen", "bishop1", "bishop2", "knight1", "knight2", "rook1", "rook2"]
N = len(PIECES)  # 8 pieces

def cell(r, c):
    return r * COLS + c

def pos_to_rc(pos):
    return pos // COLS, pos % COLS

def color(r, c):
    """Square color: 0=light, 1=dark"""
    return (r + c) % 2

def attacks_king(r1, c1, r2, c2):
    return max(abs(r1 - r2), abs(c1 - c2)) == 1

def attacks_queen(r1, c1, r2, c2):
    return (r1 == r2 or c1 == c2 or
            abs(r1 - r2) == abs(c1 - c2))

def attacks_bishop(r1, c1, r2, c2):
    return abs(r1 - r2) == abs(c1 - c2)

def attacks_rook(r1, c1, r2, c2):
    return r1 == r2 or c1 == c2

def attacks_knight(r1, c1, r2, c2):
    dr, dc = abs(r1 - r2), abs(c1 - c2)
    return (dr == 1 and dc == 2) or (dr == 2 and dc == 1)

# Map piece index to its attack function
ATTACK_FN = {
    0: attacks_king,
    1: attacks_queen,
    2: attacks_bishop,
    3: attacks_bishop,
    4: attacks_knight,
    5: attacks_knight,
    6: attacks_rook,
    7: attacks_rook,
}

def solve_all():
    solutions = []

    # Each piece gets an integer position in [0, TOTAL_CELLS)
    pos = [Int(f"pos_{p}") for p in PIECES]

    s = Solver()

    # All positions in valid range
    for p in pos:
        s.add(p >= 0, p < TOTAL_CELLS)

    # All pieces on distinct squares
    s.add(Distinct(*pos))

    # Bishops on opposite colors
    # bishop1 = pos[2], bishop2 = pos[3]
    # color of pos = ((pos // COLS) + (pos % COLS)) % 2
    b1_color = (pos[2] / COLS + pos[2] % COLS) % 2  # Z3 integer arithmetic
    b2_color = (pos[3] / COLS + pos[3] % COLS) % 2

    # Z3 doesn't do Python //; use integer division carefully
    # We'll encode opposite colors: sum of colors == 1
    # color(pos) = (row + col) % 2 = (pos//COLS + pos%COLS) % 2
    # For Z3 int vars, use: (pos / COLS + pos % COLS) % 2
    s.add((pos[2] / COLS + pos[2] % COLS) % 2 != (pos[3] / COLS + pos[3] % COLS) % 2)

    # No piece attacks another
    # We enumerate all (i,j) pairs and for each possible (pos_i, pos_j) assignment,
    # assert that those pieces don't attack each other.
    # Since positions are Z3 vars, we use a constraint table approach:
    # For each pair (i,j), add constraints over all concrete (r1,c1,r2,c2) combos.

    # Pre-build attack lookup tables for each piece-type pair
    # This avoids slow Z3 symbolic reasoning over large disjunctions.

    # We'll encode: for all concrete assignments, if pos[i]==cell(r1,c1) and pos[j]==cell(r2,c2),
    # then piece i must not attack piece j.
    # Equivalently: NOT (pos[i]==cell(r1,c1) AND pos[j]==cell(r2,c2)) for all attacking pairs.

    # For efficiency, we build a set of "forbidden (cell_i, cell_j)" pairs per piece pair.

    for i, j in itertools.combinations(range(N), 2):
        fn_i = ATTACK_FN[i]
        fn_j = ATTACK_FN[j]
        forbidden = []
        for r1, c1 in itertools.product(range(ROWS), range(COLS)):
            for r2, c2 in itertools.product(range(ROWS), range(COLS)):
                if (r1, c1) == (r2, c2):
                    continue
                # i attacks j or j attacks i
                if fn_i(r1, c1, r2, c2) or fn_j(r2, c2, r1, c1):
                    forbidden.append((cell(r1, c1), cell(r2, c2)))

        if forbidden:
            # Add: NOT (pos[i] == ci AND pos[j] == cj) for each forbidden pair
            s.add(And([
                Not(And(pos[i] == ci, pos[j] == cj))
                for ci, cj in forbidden
            ]))

    # Enumerate all solutions
    solution_count = 0
    while s.check() == sat:
        m = s.model()
        sol = [m[pos[k]].as_long() for k in range(N)]
        solutions.append(sol)
        solution_count += 1

        # Block this solution (and symmetric duplicates via blocking exact assignment)
        s.add(Or([pos[k] != sol[k] for k in range(N)]))

        if solution_count % 10 == 0:
            print(f"  Found {solution_count} solutions so far...")

    return solutions

def print_board(sol):
    board = [["." for _ in range(COLS)] for _ in range(ROWS)]
    symbols = ["K", "Q", "B", "B", "N", "N", "R", "R"]
    for k, p in enumerate(sol):
        r, c = pos_to_rc(p)
        board[r][c] = symbols[k]
    print("  " + " ".join(str(c) for c in range(COLS)))
    for r in range(ROWS):
        print(f"{r} " + " ".join(board[r]))

def canonicalize(sol):
    """
    Produce a canonical form to deduplicate solutions where
    bishop1/bishop2, knight1/knight2, rook1/rook2 are interchangeable.
    Returns a frozenset-based key.
    """
    king_pos   = sol[0]
    queen_pos  = sol[1]
    bishops    = tuple(sorted([sol[2], sol[3]]))
    knights    = tuple(sorted([sol[4], sol[5]]))
    rooks      = tuple(sorted([sol[6], sol[7]]))
    return (king_pos, queen_pos, bishops, knights, rooks)

def main():
    print("Solving chess puzzle on 6×5 board...")
    print("Pieces: King, Queen, 2 Bishops (opposite colors), 2 Knights, 2 Rooks")
    print("Constraint: No piece attacks another\n")

    solutions = solve_all()

    # Deduplicate by canonical form
    seen = {}
    for sol in solutions:
        key = canonicalize(sol)
        if key not in seen:
            seen[key] = sol

    unique = list(seen.values())
    print(f"\nTotal raw solutions: {len(solutions)}")
    print(f"Unique solutions (deduplicating piece-pair swaps): {len(unique)}\n")

    for idx, sol in enumerate(unique, 1):
        print(f"── Solution {idx} ──")
        print_board(sol)
        labels = ["King", "Queen", "Bishop1", "Bishop2", "Knight1", "Knight2", "Rook1", "Rook2"]
        for k, p in enumerate(sol):
            r, c = pos_to_rc(p)
            col_name = "light" if color(r, c) == 0 else "dark"
            sq = f"({r},{c})"
            extra = f" [{col_name}]" if "Bishop" in labels[k] else ""
            print(f"  {labels[k]}: {sq}{extra}")
        print()

if __name__ == "__main__":
    main()

Related post: Lessons Learned With the Z3 SAT/SMT Solver

Formalizing a ring theorem with Lean 4 and Claude

I’ve been testing Claude’s ability to generate Lean 4 code to prove theorems. I’ve written about a couple experiments that verified calculations. I did not write about my failed attempt to get Claude to formalize a proof of the pqr theorem for seminorms.

This time I asked Claude to formally prove the theorem from the previous post on partial fraction decomposition in an abstract setting. Unlike the previous experiments, I did not give it a proof or calculation to verify; I asked it to generate a proof.

This was the prompt:

Write a proof in Lean 4 for the following generalization of partial fraction decomposition.

If a ring R is a principle ideal domain (PID), then every element c of the field of fractions over R can be written in the form \[c = \sum_i \frac{a_i}{p_i^{r_i}}\]

where the p_i are nonassociate irreducible elements of R, the r_i are non-negative integers, and the elements a_i and p_i are relatively prime.

This was Claude’s immediate response:

I’ll work through this carefully, as it’s a substantial formalization task. Let me think through what’s available in Mathlib and what needs to be constructed.

This is a genuinely hard Lean 4 formalization problem. Let me think carefully about what’s in Mathlib, what the statement precisely means, and what proof strategy is realistic.

It also stated its overall proof strategy:

Proof strategy: Write c = a/b in lowest terms. Factor b = ∏ pᵢ^rᵢ. Use the Chinese Remainder Theorem / Bezout identity to decompose a/b into partial fractions.

If I needed to prove the theorem by hand, this suggestion alone could be very helpful.

The first attempt generated code that would not build. As before, a lot of the problems seem to stem from things moving around in Mathlib due to refactorings. There are also occasional hallucinations.

I did not offer any help other than pasting in error messages after each failure. When working on a real problem, I’d be more interactive. But as a test, I just passively pasted error messages and tried again. It took 11 iterations, but I finally got a proof. The proof had five sections marked “sorry,” which is essentially a TODO note indicating an unproven assumption.

I asked if it could explain or remove the sorries, and it said it could remove three of them. But the generated code did not build and I hit a session limit before I could debug the code. Here is the code with five sorry sections.

/-
  PartialFractions.lean  (v11 — final)
-/

import Mathlib

open UniqueFactorizationMonoid

variable {R : Type*} [CommRing R] [IsDomain R] [IsPrincipalIdealRing R]

structure PFDSummand (R : Type*) [CommRing R] where
  prime    : R
  exponent : ℕ
  numer    : R

noncomputable def PFDSummand.val [CommRing R] [IsDomain R] (s : PFDSummand R) :
    FractionRing R :=
  algebraMap R (FractionRing R) s.numer /
  algebraMap R (FractionRing R) (s.prime ^ s.exponent)

structure IsPFD [CommRing R] [IsDomain R] [IsPrincipalIdealRing R]
    (terms : List (PFDSummand R)) (c : FractionRing R) : Prop where
  irred    : ∀ s ∈ terms, Irreducible s.prime
  nonAssoc : ∀ s ∈ terms, ∀ t ∈ terms, s ≠ t → ¬Associated s.prime t.prime
  coprime  : ∀ s ∈ terms, IsCoprime s.numer s.prime
  sum_eq   : (terms.map PFDSummand.val).sum = c

private lemma algMap_ne_zero {a : R} (h : a ≠ 0) :
    (algebraMap R (FractionRing R)) a ≠ 0 :=
  map_ne_zero_iff _ (IsFractionRing.injective R (FractionRing R)) |>.mpr h

private lemma exists_num_denom (c : FractionRing R) :
    ∃ a b : R, b ≠ 0 ∧
      (algebraMap R (FractionRing R)) a / (algebraMap R (FractionRing R)) b = c := by
  obtain ⟨a, b, hb, hab⟩ := IsFractionRing.div_surjective (A := R) c
  exact ⟨a, b, nonZeroDivisors.ne_zero hb, hab⟩

-- ======================================================================
-- §6  Bézout splitting
-- ======================================================================

private lemma bezout_split (x y : R) (hcop : IsCoprime x y)
    (hx : x ≠ 0) (hy : y ≠ 0) (e : R) :
    ∃ s t : R, s * x + t * y = 1 ∧
      (algebraMap R (FractionRing R) e) /
        ((algebraMap R (FractionRing R) x) * (algebraMap R (FractionRing R) y)) =
      (algebraMap R (FractionRing R) (e * t)) / (algebraMap R (FractionRing R) x) +
      (algebraMap R (FractionRing R) (e * s)) / (algebraMap R (FractionRing R) y) := by
  obtain ⟨s, t, hst⟩ := hcop
  refine ⟨s, t, hst, ?_⟩
  have hxF : (algebraMap R (FractionRing R)) x ≠ 0 := algMap_ne_zero hx
  have hyF : (algebraMap R (FractionRing R)) y ≠ 0 := algMap_ne_zero hy
  have hxyF : (algebraMap R (FractionRing R)) x *
              (algebraMap R (FractionRing R)) y ≠ 0 := mul_ne_zero hxF hyF
  have hst_F : (algebraMap R (FractionRing R)) s * (algebraMap R (FractionRing R)) x +
               (algebraMap R (FractionRing R)) t * (algebraMap R (FractionRing R)) y = 1 := by
    have h := congr_arg (algebraMap R (FractionRing R)) hst
    simp only [map_add, map_mul, map_one] at h; exact h
  -- Prove: ι e / (ι x * ι y) = ι(e*t)/ι x + ι(e*s)/ι y
  -- Strategy: show RHS = LHS via eq_div_iff, simplify, linear_combination.
  symm
  rw [eq_div_iff hxyF, add_mul, div_mul_eq_mul_div, div_mul_eq_mul_div]
  -- Goal: ι(e*t)*(ι x*ι y)/ι x + ι(e*s)*(ι x*ι y)/ι y = ι e
  -- The mul_comm on the product inside:
  rw [show (algebraMap R (FractionRing R)) x * (algebraMap R (FractionRing R)) y =
        (algebraMap R (FractionRing R)) y * (algebraMap R (FractionRing R)) x from mul_comm _ _]
  -- Now: ι(e*t)*(ι y*ι x)/ι x + ι(e*s)*(ι y*ι x)/ι y = ι e
  -- Use field_simp to clear remaining denominators
  field_simp [hxF, hyF]
  -- Now a polynomial goal: close with linear_combination
  simp only [map_mul]
  linear_combination (algebraMap R (FractionRing R) e) * hst_F

-- ======================================================================
-- §7  Coprime numerator reduction
-- ======================================================================

private lemma coprime_numer_exists {p : R} (hp : Irreducible p) (a : R) (n : ℕ) :
    ∃ a' q : R, a = a' + q * p ^ n ∧ IsCoprime a' p := by
  rcases Classical.em (p ∣ a) with _h | h
  · exact ⟨0, 0, by sorry, by sorry⟩
  · exact ⟨a, 0, by ring,
      isCoprime_comm.mpr ((irreducible_iff_prime.mp hp).coprime_iff_not_dvd.mpr h)⟩

-- ======================================================================
-- §8  Main theorem
-- ======================================================================

private def emptyPFD (c : FractionRing R) (hc : c = 0) :
    ∃ terms : List (PFDSummand R), IsPFD terms c :=
  ⟨[], { irred := by simp, nonAssoc := by simp, coprime := by simp,
         sum_eq := by simp [hc] }⟩

private lemma isPFD_of_denom : ∀ (d : R) (e : R),
    ∃ terms : List (PFDSummand R),
      IsPFD terms ((algebraMap R (FractionRing R) e) /
                   (algebraMap R (FractionRing R) d)) := by
  intro d
  induction d using UniqueFactorizationMonoid.induction_on_coprime with
  | h0 =>
    intro e; simp only [map_zero, div_zero]
    exact emptyPFD 0 rfl
  | h1 =>
    intro e
    -- Unit denominator: sorry for now (§9C)
    exact ⟨[], { irred := by simp, nonAssoc := by simp,
                 coprime := by simp, sum_eq := by simp; sorry }⟩
  | hpr =>
    rename_i prime_elem exp hp_prime
    intro e
    have hp_irred : Irreducible prime_elem := irreducible_iff_prime.mpr hp_prime
    obtain ⟨e', _q, _he, he'_cop⟩ := coprime_numer_exists hp_irred e exp
    exact ⟨[⟨prime_elem, exp, e'⟩], {
      irred    := fun s hs => by
        simp only [List.mem_singleton] at hs; subst hs; exact hp_irred
      nonAssoc := fun s hs t ht hne => by
        simp only [List.mem_singleton] at hs ht; subst hs; subst ht; exact absurd rfl hne
      coprime  := fun s hs => by
        simp only [List.mem_singleton] at hs; subst hs; exact he'_cop
      sum_eq   := by
        simp only [List.map_singleton, List.sum_singleton, PFDSummand.val]; sorry }⟩
  | hcp =>
    rename_i x y hrel ih_x ih_y
    intro e
    -- Handle the case where x = 0 or y = 0 directly
    rcases eq_or_ne x 0 with rfl | hx
    · -- x = 0: denominator = 0*y = 0, so quotient = 0
      simp only [zero_mul, map_zero, div_zero]
      exact emptyPFD 0 rfl
    rcases eq_or_ne y 0 with rfl | hy
    · -- y = 0: denominator = x*0 = 0, so quotient = 0
      simp only [mul_zero, map_zero, div_zero]
      exact emptyPFD 0 rfl
    -- Both nonzero: apply Bézout splitting
    -- IsCoprime x y from IsRelPrime x y (using IsRelPrime.isCoprime)
    have hcop : IsCoprime x y := hrel.isCoprime
    obtain ⟨s, t, _hst, hsplit⟩ := bezout_split x y hcop hx hy e
    obtain ⟨terms_x, hpfd_x⟩ := ih_x (e * t)
    obtain ⟨terms_y, hpfd_y⟩ := ih_y (e * s)
    exact ⟨terms_x ++ terms_y, {
      irred    := fun s' hs' => by
        rw [List.mem_append] at hs'
        exact hs'.elim (hpfd_x.irred s') (hpfd_y.irred s')
      nonAssoc := fun s' hs' t' ht' hne => by
        rw [List.mem_append] at hs' ht'
        rcases hs' with hs' | hs' <;> rcases ht' with ht' | ht'
        · exact hpfd_x.nonAssoc s' hs' t' ht' hne
        · intro; sorry
        · intro; sorry
        · exact hpfd_y.nonAssoc s' hs' t' ht' hne
      coprime  := fun s' hs' => by
        rw [List.mem_append] at hs'
        exact hs'.elim (hpfd_x.coprime s') (hpfd_y.coprime s')
      sum_eq   := by
        rw [List.map_append, List.sum_append, hpfd_x.sum_eq, hpfd_y.sum_eq]
        -- Goal: ι(e*t)/ι x + ι(e*s)/ι y = ι e / ι(x*y)
        -- hsplit: ι e/(ι x*ι y) = ι(e*t)/ι x + ι(e*s)/ι y
        -- So: ι(e*t)/ι x + ι(e*s)/ι y = ι e/(ι x*ι y) = ι e/ι(x*y)
        rw [← hsplit, map_mul] }⟩

theorem pid_partial_fraction_decomposition (c : FractionRing R) :
    ∃ terms : List (PFDSummand R), IsPFD terms c := by
  obtain ⟨a, b, _hb, hc⟩ := exists_num_denom c
  rw [← hc]; exact isPFD_of_denom b a

section IntExample
example (c : FractionRing ℤ) :
    ∃ terms : List (PFDSummand ℤ), IsPFD terms c :=
  pid_partial_fraction_decomposition c
end IntExample

Partial fraction decomposition

Nearly everyone who has seen partial fraction decomposition was introduced to it as a way to compute integrals. If P(x) and Q(x) are polynomials, then you can break their ratio P(x)/Q(x) into a sum of terms that can each be integrated in closed form. As with most topics in a calculus class, partial fractions go by in a blur.

This post will look at partial fractions more generally.

Computation

Every polynomial with real coefficients can be factored into a product of linear and irreducible quadratic terms. But actually calculating this factorization is difficult if the degree of the denominator is large.

The quadratic equation is easy to use. There are analogs for 3rd and 4th order polynomials, but they’re cumbersome. And there is no formula in general for finding roots of polynomials of degree 5 or higher.

You could find the roots numerically, but if you’re going to go that route, maybe you should evaluate your integral numerically.

Still, it is useful in proving theorems to know that a partial fraction decomposition exists, even if in practice you cannot calculate it.

Complex numbers

Rational polynomials over the real numbers can be factored into powers of linear terms and irreducible quadratic terms. There are no irreducible quadratics over the complex numbers thanks to the Fundamental Theorem of Algebra, and every polynomial can be factored into a product of linear terms.

This means every rational in z can be broken into a sum of a polynomial in z and polynomials in 1/(zzi) where the zi are the roots of the denominator. This fact is important, for example, in contour integration.

Principle ideal domains

The concept of partial fraction decomposition can be generalized to the field of fractions over a ring R [1].

If the ring R is a principle ideal domain (PID) [2], then every element c of the field K of fractions over R can be written in the form

c = \sum_i \frac{a_i}{p_i^{r_i}}

where the pi are nonassociate [3] irreducible elements of R, the ri are non-negative integers, and the elements ai and pi are relatively prime.

When R is the ring of of polynomials over a field, R is a PID, and the field of fractions is the set of rational functions over that field. When the field is the real or complex numbers, we get the results above. But the field could be something else, such as a finite field.

Integers

When R is the ring of integers, the irreducible elements are prime numbers. The nonassociate condition means you can’t count p and −p as distinct elements, so practically this means we only look at positive primes. The field of fractions is the rational numbers. So the theorem above says that every rational number can be written as a sum of fractions where the denominators of the fractions are prime powers and the numerators are relatively prime to the denominators.

The way you would decompose a rational number into fractions with prime power denominators is analogous to the way you’d do partial fraction decomposition in a calculus class. For example, suppose we want to decompose 46/75. The distinct prime factors of 75 are 3 and 5, and so we’d look for fractions with denominators 3, 5, and 25, and in fact

\frac{46}{75} = \frac{1}{3} + \frac{2}{5} - \frac{3}{25}

Footnotes

[1] The field of fractions over R is the set of formal terms a/b where a and b are in R and b ≠ 0. Operations are defined by analogy with rational numbers. If R is an integral domain, the field of fractions really is a field.

[2] A ring is a PID if every ideal can be generated by a single element.

[3] Two elements of an integral domain are said to be associate if they generate the same ideal.