//
sign in
Profile
by @danabra.mov
Profile
by @dansshadow.bsky.social
Profile
by @jimpick.com
AviHandle
by @danabra.mov
AviHandle
by @dansshadow.bsky.social
AviHandle
by @katherine.computer
EventsList
by @katherine.computer
ProfileHeader
by @dansshadow.bsky.social
ProfileHeader
by @danabra.mov
ProfileMedia
by @danabra.mov
ProfilePlays
by @danabra.mov
ProfilePosts
by @danabra.mov
ProfilePosts
by @dansshadow.bsky.social
ProfileReplies
by @danabra.mov
Record
by @atsui.org
Skircle
by @danabra.mov
StreamPlacePlaylist
by @katherine.computer
+ new component
ProfilePosts









Loading...
Placing all six chess pieces on a 6 by 5 board so that no piece attacks another. Solution using Claude to generate Z3/Python. https://www.johndcook.com/blog/2026/06/20/z3-python-claude/
Pentagonal numbers and drawing inspiration from bug fixes https://johndcookconsulting.substack.com/p/pentagonal-numbers
Formalizing a ring theorem with Lean 4 and Claude https://www.johndcook.com/blog/2026/06/17/rings-with-lean-claude/
Pulling on a thread: How an approximation lead to a series of nine posts. https://www.johndcook.com/blog/2026/06/10/pulling-on-a-thread/
Quaternion Rotations, Claude, and Lean https://www.johndcook.com/blog/2026/06/15/quaternions-claude-lean/
21h
5d
4d
Partial fraction decomposition for more than calculating integrals https://www.johndcook.com/blog/2026/06/16/partial-fraction-decomposition/
Solving a chess puzzle with Claude and Prolog https://www.johndcook.com/blog/2026/06/11/prolog-claude/
Testing pentagonal numbers https://www.johndcook.com/blog/2026/06/15/testing-pentagonal-numbers/
Formally proving a calculation with Claude and Lean https://www.johndcook.com/blog/2026/06/10/claude-and-lean/
Solving another chess puzzle using AI-generated Prolog code https://www.johndcook.com/blog/2026/06/15/writing-prolog-with-chatgpt/
11d
5d
10d
5d
5d
10d
5d
John D. Cook
John D. Cook
John D. Cook
John D. Cook
John D. Cook
John D. Cook
John D. Cook
John D. Cook
John D. Cook
John D. Cook
Inspiration from bug fixes
johndcookconsulting.substack.com
Pentagonal Numbers
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
www.johndcook.com
I ran an experiment today to see whether Claude could generate LEAN code to prove a calculation at the bottom of this post, six lines of calculus. I started with this prompt This page contains a mathematical proof that a Fourier coefficient, a_n, is given in terms of a Bessel function. The LaTeX source for the SVG image is contained in the alt tag of the image. Generate a formal proof of the result using LEAN. and give it the URL of the post. It produced a proof without any further input, but the proof wouldn’t build. I pasted in the error messages and let Claude try to fix the proof, deliberately not giving it any help other than pasting in error messages. It took eight iterations, but eventually it worked. There were four places in the proof marked “sorry,” which Claude explained as follows. **Why the`sorry`s?** Real.besselJ definition or a `Kepler` namespace. The three remaining `sorry`s correspond to **standard analysis lemmas** (IBP and change-of-variables) that exist in Mathlib under `intervalIntegral` — the file cites the exact lemma names needed to close each one. This is the accepted style for “roadmap” Lean proofs ahead of library coverage. The full proof is given below. /- KeplerBessel.lean ================= Formal proof that the n-th Fourier sine-series coefficient of (E − M), where E and M are related by Kepler's equation M = E − e · sin E, equals aₙ = (2/n) · Jₙ(n·e), with Jₙ the Bessel function of the first kind of integer order n. Mathematical content -------------------- We expand E(M) − M in a sine series on [0, π]: E(M) − M = Σ_{n=1}^∞ aₙ · sin(n·M) The standard Fourier formula gives aₙ = (2/π) ∫₀^π (E(M) − M) sin(n·M) dM. Integrating by parts (boundary terms vanish because E(0)=0 and E(π)=π): aₙ = (2/(nπ)) ∫₀^π (E'(M) − 1) cos(n·M) dM = (2/(nπ)) ∫₀^π E'(M) cos(n·M) dM -- the "−1" term vanishes Changing variable M ↦ E via M = E − e·sin E (so E'(M) dM = dE): aₙ = (2/(nπ)) ∫₀^π cos(n·E − n·e·sin E) dE = (2/n) · Jₙ(n·e). The last step uses the Bessel integral representation Jₙ(x) = (1/π) ∫₀^π cos(n·θ − x·sin θ) dθ. -/ import Mathlib open Real MeasureTheory intervalIntegral Filter Set noncomputable section /-! --------------------------------------------------------------- §1 Variables --------------------------------------------------------------- -/ variable (e : ℝ) (he : 0 ≤ e) (he1 : e < 1) /-! --------------------------------------------------------------- §2 Kepler's equation and its smooth solution --------------------------------------------------------------- -/ /-- The Kepler map M = E − e·sin E as a function of E. -/ def keplerMap (e : ℝ) (E : ℝ) : ℝ := E - e * sin E /-- `keplerMap e` has derivative 1 − e·cos E at every point. -/ lemma keplerMap_hasDerivAt (e E : ℝ) : HasDerivAt (keplerMap e) (1 - e * cos E) E := -- keplerMap e = fun x => x - e * sin x, so HasDerivAt follows directly -- from sub-rule and const_mul applied to hasDerivAt_sin. (hasDerivAt_id E).sub ((hasDerivAt_sin E).const_mul e) /-- The derivative of `keplerMap e` is positive when e < 1. -/ lemma keplerMap_deriv_pos {e' : ℝ} (he' : 0 ≤ e') (he1' : e' < 1) (E : ℝ) : 0 < 1 - e' * cos E := by have hcos : cos E ≤ 1 := cos_le_one E nlinarith [mul_le_of_le_one_right he' hcos] /-- `keplerMap e` is strictly monotone when e < 1. Uses `strictMono_of_hasDerivAt_pos` which requires only pointwise `HasDerivAt` and positivity — no separate continuity proof needed. -/ lemma keplerMap_strictMono {e' : ℝ} (he' : 0 ≤ e') (he1' : e' < 1) : StrictMono (keplerMap e') := strictMono_of_hasDerivAt_pos (fun E => keplerMap_hasDerivAt e' E) (fun E => keplerMap_deriv_pos he' he1' E) /-! We axiomatise the inverse eccAnom : ℝ → ℝ → ℝ and its key properties, all of which follow from the Inverse Function Theorem applied to the smooth, strictly monotone map keplerMap e. -/ /-- The eccentric anomaly: the smooth right-inverse of `keplerMap e`. -/ axiom eccAnom (e : ℝ) : ℝ → ℝ /-- `eccAnom e M` satisfies Kepler's equation. -/ axiom eccAnom_kepler (e M : ℝ) : keplerMap e (eccAnom e M) = M /-- `eccAnom e` is differentiable, derivative = 1/(1 − e·cos(eccAnom e M)). -/ axiom eccAnom_hasDerivAt (e M : ℝ) : HasDerivAt (eccAnom e) (1 / (1 - e * cos (eccAnom e M))) M /-- Boundary value at 0. -/ axiom eccAnom_zero (e : ℝ) : eccAnom e 0 = 0 /-- Boundary value at π. -/ axiom eccAnom_pi (e : ℝ) : eccAnom e π = π /-! --------------------------------------------------------------- §3 Bessel function of the first kind (integer order) Defined by the classical integral representation. --------------------------------------------------------------- -/ /-- Bessel function J_n(x) via its integral representation. -/ def besselJ (n : ℕ) (x : ℝ) : ℝ := (1 / π) * ∫ θ in (0 : ℝ)..π, cos (↑n * θ - x * sin θ) /-! --------------------------------------------------------------- §4 Fourier coefficient Named keplerFourierCoeff to avoid clashing with Mathlib's own `fourierCoeff` which is defined on AddCircle. --------------------------------------------------------------- -/ /-- The n-th Fourier sine coefficient of eccAnom e M − M on [0,π]. -/ def keplerFourierCoeff (e : ℝ) (n : ℕ) : ℝ := (2 / π) * ∫ M in (0 : ℝ)..π, (eccAnom e M - M) * sin (↑n * M) /-! --------------------------------------------------------------- §5 Main theorem --------------------------------------------------------------- -/ /-- **Main theorem.** For n ≥ 1, the Fourier sine coefficient of the eccentric-anomaly displacement satisfies aₙ = (2/n) · Jₙ(n·e). -/ theorem keplerFourierCoeff_eq_besselJ (n : ℕ) (hn : 1 ≤ n) : keplerFourierCoeff e n = (2 / (n : ℝ)) * besselJ n (↑n * e) := by simp only [keplerFourierCoeff, besselJ] -- Goal: -- (2/π) · ∫₀^π (E(M)−M)·sin(nM) dM -- = (2/n) · (1/π) · ∫₀^π cos(nθ − ne·sinθ) dθ -- ── Step 1: Integration by parts ───────────────────────────────────── -- u = E(M)−M, dv = sin(nM)dM → v = −cos(nM)/n -- Boundary: [uv]₀^π = 0 by eccAnom_zero, eccAnom_pi. -- Result: (2/π)∫(E−M)sin(nM)dM = (2/π)(1/n)∫(E'(M)−1)cos(nM)dM -- -- Mathlib lemma: intervalIntegral.integral_mul_deriv -- (or integral_deriv_mul_eq_sub_of_hasDerivAt applied to -- u = eccAnom e − id, v = −sin(n··)/n) have step1 : (2 / π) * ∫ M in (0 : ℝ)..π, (eccAnom e M - M) * sin (↑n * M) = (2 / π) * (1 / ↑n) * ∫ M in (0 : ℝ)..π, (deriv (eccAnom e) M - 1) * cos (↑n * M) := by sorry -- ── Step 2: The "−1" integral vanishes ─────────────────────────────── -- ∫₀^π cos(nM) dM = [sin(nM)/n]₀^π = 0 (integer n ≥ 1) -- Mathlib: integral_cos, Real.sin_nat_mul_pi have cos_int_zero : ∫ M in (0 : ℝ)..π, cos (↑n * M) = 0 := by sorry have step2 : ∫ M in (0 : ℝ)..π, (deriv (eccAnom e) M - 1) * cos (↑n * M) = ∫ M in (0 : ℝ)..π, deriv (eccAnom e) M * cos (↑n * M) := by have key : ∀ M : ℝ, (deriv (eccAnom e) M - 1) * cos (↑n * M) = deriv (eccAnom e) M * cos (↑n * M) - cos (↑n * M) := by intro M; ring simp_rw [key] rw [intervalIntegral.integral_sub _ _] · rw [cos_int_zero, sub_zero] · -- IntervalIntegrable (deriv (eccAnom e) · cos(n··)) sorry · exact (continuous_cos.comp (continuous_const.mul continuous_id')).intervalIntegrable 0 π -- ── Step 3: Change of variable M ↦ E via Kepler's equation ─────────── -- Under M = E − e·sin E: E'(M) dM = dE, cos(nM) = cos(nE − ne·sinE) -- Mathlib: MeasureTheory.integral_image_eq_integral_abs_deriv_smul -- or intervalIntegral.integral_comp_deriv have step3 : ∫ M in (0 : ℝ)..π, deriv (eccAnom e) M * cos (↑n * M) = ∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E) := by sorry -- ── Step 4: Recognise the Bessel integral ──────────────────────────── -- (1/π)∫₀^π cos(nE − ne·sinE)dE = J_n(ne) by definition. have step4 : (1 / π) * ∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E) = besselJ n (↑n * e) := by simp only [besselJ] -- ── Combine all steps ──────────────────────────────────────────────── have hpi : π > 0 := Real.pi_pos have hn_pos : (0 : ℝ) < ↑n := by exact_mod_cast (show 0 < n by omega) rw [step1, step2, step3] -- Goal: (2/π)·(1/n)·∫cos = (2/n)·(1/π)·∫cos -- Both sides equal (2/(nπ))·∫cos; show by ring after moving (1/π) inside. have factored : (2 / π) * (1 / ↑n) * ∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E) = (2 / ↑n) * ((1 / π) * ∫ E in (0 : ℝ)..π, cos (↑n * E - ↑n * e * sin E)) := by ring rw [factored, step4] end /-! --------------------------------------------------------------- §6 Summary of proof obligations Axioms (all standard consequences of the Inverse Function Theorem applied to the smooth, strictly increasing map keplerMap e): • eccAnom — existence of the smooth inverse • eccAnom_kepler — it satisfies Kepler's equation • eccAnom_hasDerivAt — C¹ with derivative 1/(1 − e·cos E) • eccAnom_zero — E(0) = 0 • eccAnom_pi — E(π) = π `sorry`s (each closes with a standard Mathlib lemma): • step1 IBP via intervalIntegral.integral_deriv_mul_eq_sub_of_hasDerivAt • cos_int_zero ∫₀^π cos(nM)dM = 0 via integral_cos + Real.sin_nat_mul_pi • step2 integrability IntervalIntegrable for deriv(eccAnom e)·cos(n·) • step3 Change of variables via MeasureTheory.integral_image_eq_integral_abs_deriv_smul §7 Finding minimal imports Once the file builds cleanly, add at the bottom: #min_imports and the Lean Infoview will report the exact minimal import list for the version of Mathlib you have installed. --------------------------------------------------------------- -/
www.johndcook.com
Formally proving a calculation with Claude and LEAN
Nearly everyone who as 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/(_z_ − _z_ _i_) where the _z_ _i_ 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 where the _p_ _i_ are nonassociate [3] irreducible elements of _R_ , the _r_ _i_ are non-negative integers, and the elements _a_ _i_ and _p_ _i_ 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 ## 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 an 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.
www.johndcook.com
Partial fraction decomposition
Often there’s a thread running through a sequence of my posts. Sometimes I make this explicit and sometimes I don’t. The latest thread started with this post commenting on a tweet that observed that exp(− _x_ ²) ≈ (1 + cos(sin(_x_ + _x_))/2. Some people said online that that the approximation is simply due to the first few terms of the Taylor series on both sides matching up, so I wrote a follow up post explaining that it’s not that simple. The series for the left hand side alternates and converges very slowly, which lead to the post on naively summing an alternating series. The series for the right hand side lead to this point on partitions over permutations. Integrating the right hand side lead to this post on how the simplest numerical integration rule works shockingly well on some problems. The exact value of the integral turns out to be given by a Bessel function, details given in this post. Mr. Bessel’s interest in the functions now named after him started with looking closely at a solution to Kepler’s equation in orbital mechanics. Thinking about Kepler’s equation lead to the posts on the Laplace limit and on series acceleration. I may be done pulling on this thread. I don’t have anything else in mind that I want to explore for now, but you never know.
www.johndcook.com
Pulling on a thread