[{"id":247152,"date":"2026-06-11T08:14:26","date_gmt":"2026-06-11T13:14:26","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247152"},"modified":"2026-06-11T08:14:26","modified_gmt":"2026-06-11T13:14:26","slug":"prolog-claude","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/06\/11\/prolog-claude\/","title":{"rendered":"Solving a chess puzzle with Claude and Prolog"},"content":{"rendered":"

Prolog is the original logic programming language. The name comes from pro<\/span>gramming in log<\/span>ic. More specifically, the name comes from pro<\/span>grammation en log<\/span>ique<\/span><\/em> because the inventor of the language, Philippe Roussel, is French.<\/p>\n

Prolog has its advantages and disadvantages. One of the advantages is that the language represents logical problems directly. One of the disadvantages is that the syntax can be quirky. But if an LLM is writing the code, or at least helping to write the code, the syntax doesn’t matter so much.<\/p>\n

I wanted to see how well Claude (Sonnet 4.6, medium effort) could solve a chess puzzle<\/a> by Martin Gardner that I wrote about a little over a year ago. I chose a relatively obscure problem rather than something like the Eight Queens puzzle because an LLM could simply quote one of countless articles on the puzzle.<\/p>\n

The puzzle<\/h2>\n

As I stated in the post last year, the task is to place two rooks, two bishops, and two knights on a 4 by 4 chessboard so that no piece attacks any other.<\/p>\n

<\/p>\n

There are two basic solutions, twelve if you count rotations and reflections as different solutions.<\/p>\n

Prolog results<\/h2>\n

Claude wrote an SWI-Prolog program that I ran with<\/p>\n

\r\nswipl -g \"run, halt\" chess_placement.pl\r\n<\/pre>\n
and it gave the following output.<\/p>\n
\r\nTotal solutions (raw): 12\r\nSolutions up to D4 symmetry: 2\r\nSolutions up to D4 symmetry: 2\r\n\r\nSolution 1:\r\nB . . N\r\nN . . B\r\n. R . .\r\n. . R .\r\n\r\nSolution 2:\r\n. B N .\r\n. . . R\r\nR . . .\r\n. N B .\r\n<\/pre>\n
Claude made a nice image of the solutions.<\/p>\n
<\/p>\n
Note that you can make eight variations of the first solution by a combination of rotations and flips. There are only four variations of the second solution because two rotations put you back where you started.<\/p>\n
Prolog source<\/h2>\n
Here’s the Prolog code that Claude generated.<\/p>\n
\r\n% ============================================================\r\n%  Non-attacking placement of 2 Rooks, 2 Bishops, 2 Knights\r\n%  on a 4x4 chessboard\r\n%\r\n%  Run with SWI-Prolog:\r\n%    swipl -g \"run, halt\" chess_placement.pl\r\n%\r\n%  Or interactively:\r\n%    swipl chess_placement.pl\r\n%    ?- run.\r\n%    ?- count_solutions(N).\r\n%    ?- count_unique(N).\r\n% ============================================================\r\n\r\n:- use_module(library(lists)).\r\n\r\n% ---- Board squares: (Row, Col) with Row, Col in 0..3 ----------------\r\n\r\nsquare(R, C) :- member(R, [0,1,2,3]), member(C, [0,1,2,3]).\r\n\r\n% ---- Attack predicates ----------------------------------------------\r\n\r\nrook_attacks(R,_,R,_).          % same row\r\nrook_attacks(_,C,_,C).          % same column\r\n\r\nbishop_attacks(R1,C1,R2,C2) :-\r\n    DR is abs(R1-R2), DC is abs(C1-C2),\r\n    DR =:= DC, DR > 0.\r\n\r\nknight_attacks(R1,C1,R2,C2) :-\r\n    DR is abs(R1-R2), DC is abs(C1-C2),\r\n    ( DR =:= 1, DC =:= 2 ; DR =:= 2, DC =:= 1 ).\r\n\r\n% attacks(+Type1, +R1, +C1, +Type2, +R2, +C2)\r\n% True if either piece threatens the other's square.\r\nattacks(rook,   R1,C1, _,     R2,C2) :- rook_attacks(R1,C1,R2,C2).\r\nattacks(_,      R1,C1, rook,  R2,C2) :- rook_attacks(R1,C1,R2,C2).\r\nattacks(bishop, R1,C1, _,     R2,C2) :- bishop_attacks(R1,C1,R2,C2).\r\nattacks(_,      R1,C1, bishop,R2,C2) :- bishop_attacks(R1,C1,R2,C2).\r\nattacks(knight, R1,C1, _,     R2,C2) :- knight_attacks(R1,C1,R2,C2).\r\nattacks(_,      R1,C1, knight,R2,C2) :- knight_attacks(R1,C1,R2,C2).\r\n\r\n% ---- All pairs safe -------------------------------------------------\r\n\r\n% Check all unordered pairs in a list.\r\nno_attack_pair(_, []).\r\nno_attack_pair(T1-(R1,C1), [T2-(R2,C2)|Rest]) :-\r\n    \\+ attacks(T1,R1,C1,T2,R2,C2),\r\n    no_attack_pair(T1-(R1,C1), Rest).\r\n\r\nall_pairs_safe([]).\r\nall_pairs_safe([P|Rest]) :-\r\n    no_attack_pair(P, Rest),\r\n    all_pairs_safe(Rest).\r\n\r\n% ---- Generate a placement -------------------------------------------\r\n%\r\n%  Placement = [rook-(R1,C1), rook-(R2,C2),\r\n%               bishop-(R3,C3), bishop-(R4,C4),\r\n%               knight-(R5,C5), knight-(R6,C6)]\r\n%\r\n%  We enforce an ordering within each pair to avoid counting\r\n%  same-type duplicates (rook1 vs rook2 swapped).\r\n%  Squares are represented canonically: (R,C) with R*4+C as index.\r\n\r\nsq_index(R, C, I) :- I is R*4 + C.\r\n\r\nplacement(Placement) :-\r\n    % --- Rooks ---\r\n    square(R1r, C1r), square(R2r, C2r),\r\n    sq_index(R1r,C1r,I1r), sq_index(R2r,C2r,I2r),\r\n    I1r < I2r,          % canonical order for the two rooks\r\n\r\n    % --- Bishops ---\r\n    square(R1b, C1b), square(R2b, C2b),\r\n    sq_index(R1b,C1b,I1b), sq_index(R2b,C2b,I2b),\r\n    I1b < I2b,          % canonical order for the two bishops\r\n\r\n    % --- Knights ---\r\n    square(R1n, C1n), square(R2n, C2n),\r\n    sq_index(R1n,C1n,I1n), sq_index(R2n,C2n,I2n),\r\n    I1n < I2n,          % canonical order for the two knights\r\n\r\n    % --- All six squares distinct ---\r\n    Squares = [(R1r,C1r),(R2r,C2r),(R1b,C1b),(R2b,C2b),(R1n,C1n),(R2n,C2n)],\r\n    all_distinct_squares(Squares),\r\n\r\n    % --- Build placement list ---\r\n    Placement = [ rook-(R1r,C1r),   rook-(R2r,C2r),\r\n                  bishop-(R1b,C1b), bishop-(R2b,C2b),\r\n                  knight-(R1n,C1n), knight-(R2n,C2n) ],\r\n\r\n    % --- No piece attacks any other ---\r\n    all_pairs_safe(Placement).\r\n\r\nall_distinct_squares([]).\r\nall_distinct_squares([S|Rest]) :-\r\n    \\+ member(S, Rest),\r\n    all_distinct_squares(Rest).\r\n\r\n% ---- Count all solutions --------------------------------------------\r\n\r\ncount_solutions(N) :-\r\n    findall(P, placement(P), Ps),\r\n    length(Ps, N),\r\n    format(\"Total solutions (raw): ~w~n\", [N]).\r\n\r\n% ---- D4 symmetry transformations on a 4x4 board ---------------------\r\n%  Squares (R,C) in 0..3\r\n\r\ntransform(identity, R, C, R,  C).\r\ntransform(rot90,    R, C, C,  Nr) :- Nr is 3-R.\r\ntransform(rot180,   R, C, Nr, Nc) :- Nr is 3-R, Nc is 3-C.\r\ntransform(rot270,   R, C, Nc, R)  :- Nc is 3-C.\r\ntransform(flipH,    R, C, R,  Nc) :- Nc is 3-C.\r\ntransform(flipV,    R, C, Nr, C)  :- Nr is 3-R.\r\ntransform(flipD1,   R, C, C,  R).\r\ntransform(flipD2,   R, C, Nr, Nc) :- Nr is 3-C, Nc is 3-R.\r\n\r\napply_transform(_, [], []).\r\napply_transform(T, [Type-(R,C)|Rest], [Type-(NR,NC)|TRest]) :-\r\n    transform(T, R, C, NR, NC),\r\n    apply_transform(T, Rest, TRest).\r\n\r\n% Canonical form: sort pieces within same-type pairs, then sort the\r\n% whole list to get a unique representative.\r\ncanonical_placement(Placement, Canonical) :-\r\n    findall(T, member(T,[identity,rot90,rot180,rot270,\r\n                          flipH,flipV,flipD1,flipD2]), Ts),\r\n    maplist(transform_and_sort(Placement), Ts, AllForms),\r\n    msort(AllForms, Sorted),\r\n    Sorted = [Canonical|_].\r\n\r\ntransform_and_sort(Placement, T, Sorted) :-\r\n    apply_transform(T, Placement, TPl),\r\n    msort(TPl, Sorted).\r\n\r\n% ---- Count solutions up to D4 symmetry ------------------------------\r\n\r\ncount_unique(N) :-\r\n    findall(P, placement(P), Ps),\r\n    maplist(canonical_placement, Ps, Canonicals),\r\n    list_to_set(Canonicals, Unique),\r\n    length(Unique, N),\r\n    format(\"Solutions up to D4 symmetry: ~w~n\", [N]).\r\n\r\n% ---- Pretty-print a board -------------------------------------------\r\n\r\nprint_board(Placement) :-\r\n    forall(member(R, [0,1,2,3]),\r\n           print_row(R, Placement)),\r\n    nl.\r\n\r\nprint_row(R, Placement) :-\r\n    forall(member(C, [0,1,2,3]),\r\n           print_cell(R, C, Placement)),\r\n    nl.\r\n\r\nprint_cell(R, C, Placement) :-\r\n    (   member(rook-(R,C),   Placement) -> write('R ')\r\n    ;   member(bishop-(R,C), Placement) -> write('B ')\r\n    ;   member(knight-(R,C), Placement) -> write('N ')\r\n    ;   write('. ')\r\n    ).\r\n\r\n% ---- Print all unique solutions -------------------------------------\r\n\r\nprint_unique_solutions :-\r\n    findall(P, placement(P), Ps),\r\n    maplist(canonical_placement, Ps, Canonicals),\r\n    list_to_set(Canonicals, Unique),\r\n    length(Unique, N),\r\n    format(\"~nSolutions up to D4 symmetry: ~w~n~n\", [N]),\r\n    forall(nth1(I, Unique, Sol),\r\n           ( format(\"Solution ~w:~n\", [I]),\r\n             print_board(Sol) )).\r\n\r\n% ---- Top-level entry point ------------------------------------------\r\n\r\nrun :-\r\n    count_solutions(Raw),\r\n    count_unique(Sym),\r\n    format(\"~n\"),\r\n    print_unique_solutions,\r\n    format(\"Summary: ~w raw solutions, ~w up to D4 symmetry.~n\",\r\n           [Raw, Sym]).\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
Prolog is the original logic programming language. The name comes from programming in logic. More specifically, the name comes from programmation en logique because the inventor of the language, Philippe Roussel, is French. Prolog has its advantages and disadvantages. One of the advantages is that the language represents logical problems directly. One of the disadvantages […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[5],"tags":[33,150,108],"class_list":["post-247152","post","type-post","status-publish","format-standard","hentry","category-computing","tag-chess","tag-logic","tag-programming"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247152","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247152"}],"version-history":[{"count":3,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247152\/revisions"}],"predecessor-version":[{"id":247157,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247152\/revisions\/247157"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247152"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247152"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247152"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247144,"date":"2026-06-10T18:11:01","date_gmt":"2026-06-10T23:11:01","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247144"},"modified":"2026-06-11T07:35:14","modified_gmt":"2026-06-11T12:35:14","slug":"claude-and-lean","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/06\/10\/claude-and-lean\/","title":{"rendered":"Formally proving a calculation with Claude and Lean"},"content":{"rendered":"
I ran an experiment today to see whether Claude [1] could generate Lean code to prove a calculation at the bottom of this post<\/a>, six lines of calculus.<\/p>\n
I started with this prompt<\/p>\n
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.<\/p>\n
and give it the URL of the post. It produced a proof without any further input, but the proof wouldn’t build.<\/p>\n
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.<\/p>\n
There were four places in the proof marked “sorry,” which Claude explained as follows.<\/p>\n
Why the sorry<\/code>s?<\/b><\/p>\n
Real.besselJ definition or a Kepler<\/code> namespace. The three remaining sorry<\/code>s correspond to standard analysis lemmas<\/strong> (IBP and change-of-variables) that exist in Mathlib under intervalIntegral<\/code> \u2014 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.<\/p>\n
The full proof is given below.<\/p>\n
\/-\r\n  KeplerBessel.lean\r\n  =================\r\n  Formal proof that the n-th Fourier sine-series coefficient of (E \u2212 M),\r\n  where E and M are related by Kepler's equation  M = E \u2212 e \u00b7 sin E,\r\n  equals  a\u2099 = (2\/n) \u00b7 J\u2099(n\u00b7e),\r\n  with J\u2099 the Bessel function of the first kind of integer order n.\r\n\r\n  Mathematical content\r\n  --------------------\r\n  We expand  E(M) \u2212 M  in a sine series on [0, \u03c0]:\r\n\r\n      E(M) \u2212 M = \u03a3_{n=1}^\u221e  a\u2099 \u00b7 sin(n\u00b7M)\r\n\r\n  The standard Fourier formula gives\r\n\r\n      a\u2099 = (2\/\u03c0) \u222b\u2080^\u03c0 (E(M) \u2212 M) sin(n\u00b7M) dM.\r\n\r\n  Integrating by parts (boundary terms vanish because E(0)=0 and E(\u03c0)=\u03c0):\r\n\r\n      a\u2099 = (2\/(n\u03c0)) \u222b\u2080^\u03c0 (E'(M) \u2212 1) cos(n\u00b7M) dM\r\n         = (2\/(n\u03c0)) \u222b\u2080^\u03c0 E'(M) cos(n\u00b7M) dM     -- the \"\u22121\" term vanishes\r\n\r\n  Changing variable M \u21a6 E via M = E \u2212 e\u00b7sin E  (so E'(M) dM = dE):\r\n\r\n      a\u2099 = (2\/(n\u03c0)) \u222b\u2080^\u03c0 cos(n\u00b7E \u2212 n\u00b7e\u00b7sin E) dE\r\n         = (2\/n) \u00b7 J\u2099(n\u00b7e).\r\n\r\n  The last step uses the Bessel integral representation\r\n      J\u2099(x) = (1\/\u03c0) \u222b\u2080^\u03c0 cos(n\u00b7\u03b8 \u2212 x\u00b7sin \u03b8) d\u03b8.\r\n-\/\r\n\r\nimport Mathlib\r\n\r\nopen Real MeasureTheory intervalIntegral Filter Set\r\n\r\nnoncomputable section\r\n\r\n\/-! ---------------------------------------------------------------\r\n    \u00a71  Variables\r\n    --------------------------------------------------------------- -\/\r\n\r\nvariable (e : \u211d) (he : 0 \u2264 e) (he1 : e < 1) \/-! --------------------------------------------------------------- \u00a72 Kepler's equation and its smooth solution --------------------------------------------------------------- -\/ \/-- The Kepler map M = E \u2212 e\u00b7sin E as a function of E. -\/ def keplerMap (e : \u211d) (E : \u211d) : \u211d := E - e * sin E \/-- `keplerMap e` has derivative 1 \u2212 e\u00b7cos E at every point. -\/ lemma keplerMap_hasDerivAt (e E : \u211d) : HasDerivAt (keplerMap e) (1 - e * cos E) E := -- keplerMap e = fun x => x - e * sin x, so HasDerivAt follows directly\r\n  -- from sub-rule and const_mul applied to hasDerivAt_sin.\r\n  (hasDerivAt_id E).sub ((hasDerivAt_sin E).const_mul e)\r\n\r\n\/-- The derivative of `keplerMap e` is positive when e < 1. -\/\r\nlemma keplerMap_deriv_pos {e' : \u211d} (he' : 0 \u2264 e') (he1' : e' < 1) (E : \u211d) :\r\n    0 < 1 - e' * cos E := by\r\n  have hcos : cos E \u2264 1 := cos_le_one E\r\n  nlinarith [mul_le_of_le_one_right he' hcos]\r\n\r\n\/-- `keplerMap e` is strictly monotone when e < 1.\r\n    Uses `strictMono_of_hasDerivAt_pos` which requires only pointwise\r\n    `HasDerivAt` and positivity \u2014 no separate continuity proof needed. -\/\r\nlemma keplerMap_strictMono {e' : \u211d} (he' : 0 \u2264 e') (he1' : e' < 1) : StrictMono (keplerMap e') := strictMono_of_hasDerivAt_pos (fun E => keplerMap_hasDerivAt e' E)\r\n    (fun E => keplerMap_deriv_pos he' he1' E)\r\n\r\n\/-!\r\n  We axiomatise the inverse  eccAnom : \u211d \u2192 \u211d \u2192 \u211d  and its key\r\n  properties, all of which follow from the Inverse Function Theorem\r\n  applied to the smooth, strictly monotone map  keplerMap e.\r\n-\/\r\n\r\n\/-- The eccentric anomaly: the smooth right-inverse of `keplerMap e`. -\/\r\naxiom eccAnom (e : \u211d) : \u211d \u2192 \u211d\r\n\r\n\/-- `eccAnom e M` satisfies Kepler's equation. -\/\r\naxiom eccAnom_kepler (e M : \u211d) :\r\n    keplerMap e (eccAnom e M) = M\r\n\r\n\/-- `eccAnom e` is differentiable, derivative = 1\/(1 \u2212 e\u00b7cos(eccAnom e M)). -\/\r\naxiom eccAnom_hasDerivAt (e M : \u211d) :\r\n    HasDerivAt (eccAnom e) (1 \/ (1 - e * cos (eccAnom e M))) M\r\n\r\n\/-- Boundary value at 0. -\/\r\naxiom eccAnom_zero (e : \u211d) : eccAnom e 0 = 0\r\n\r\n\/-- Boundary value at \u03c0. -\/\r\naxiom eccAnom_pi (e : \u211d) : eccAnom e \u03c0 = \u03c0\r\n\r\n\/-! ---------------------------------------------------------------\r\n    \u00a73  Bessel function of the first kind (integer order)\r\n\r\n    Defined by the classical integral representation.\r\n    --------------------------------------------------------------- -\/\r\n\r\n\/-- Bessel function J_n(x) via its integral representation. -\/\r\ndef besselJ (n : \u2115) (x : \u211d) : \u211d :=\r\n  (1 \/ \u03c0) * \u222b \u03b8 in (0 : \u211d)..\u03c0, cos (\u2191n * \u03b8 - x * sin \u03b8)\r\n\r\n\/-! ---------------------------------------------------------------\r\n    \u00a74  Fourier coefficient\r\n\r\n    Named  keplerFourierCoeff  to avoid clashing with Mathlib's own\r\n    `fourierCoeff` which is defined on  AddCircle.\r\n    --------------------------------------------------------------- -\/\r\n\r\n\/-- The n-th Fourier sine coefficient of  eccAnom e M \u2212 M  on [0,\u03c0]. -\/\r\ndef keplerFourierCoeff (e : \u211d) (n : \u2115) : \u211d :=\r\n  (2 \/ \u03c0) * \u222b M in (0 : \u211d)..\u03c0,\r\n    (eccAnom e M - M) * sin (\u2191n * M)\r\n\r\n\/-! ---------------------------------------------------------------\r\n    \u00a75  Main theorem\r\n    --------------------------------------------------------------- -\/\r\n\r\n\/--\r\n  **Main theorem.** For n \u2265 1, the Fourier sine coefficient of the\r\n  eccentric-anomaly displacement satisfies  a\u2099 = (2\/n) \u00b7 J\u2099(n\u00b7e).\r\n-\/\r\ntheorem keplerFourierCoeff_eq_besselJ (n : \u2115) (hn : 1 \u2264 n) :\r\n    keplerFourierCoeff e n = (2 \/ (n : \u211d)) * besselJ n (\u2191n * e) := by\r\n  simp only [keplerFourierCoeff, besselJ]\r\n  -- Goal:\r\n  --   (2\/\u03c0) \u00b7 \u222b\u2080^\u03c0 (E(M)\u2212M)\u00b7sin(nM) dM\r\n  --   = (2\/n) \u00b7 (1\/\u03c0) \u00b7 \u222b\u2080^\u03c0 cos(n\u03b8 \u2212 ne\u00b7sin\u03b8) d\u03b8\r\n\r\n  -- \u2500\u2500 Step 1: Integration by parts \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\r\n  -- u = E(M)\u2212M,  dv = sin(nM)dM  \u2192  v = \u2212cos(nM)\/n\r\n  -- Boundary: [uv]\u2080^\u03c0 = 0 by eccAnom_zero, eccAnom_pi.\r\n  -- Result: (2\/\u03c0)\u222b(E\u2212M)sin(nM)dM = (2\/\u03c0)(1\/n)\u222b(E'(M)\u22121)cos(nM)dM\r\n  --\r\n  -- Mathlib lemma: intervalIntegral.integral_mul_deriv\r\n  --   (or integral_deriv_mul_eq_sub_of_hasDerivAt applied to\r\n  --    u = eccAnom e \u2212 id,  v = \u2212sin(n\u00b7\u00b7)\/n)\r\n  have step1 :\r\n      (2 \/ \u03c0) * \u222b M in (0 : \u211d)..\u03c0, (eccAnom e M - M) * sin (\u2191n * M)\r\n      = (2 \/ \u03c0) * (1 \/ \u2191n) *\r\n          \u222b M in (0 : \u211d)..\u03c0, (deriv (eccAnom e) M - 1) * cos (\u2191n * M) := by\r\n    sorry\r\n\r\n  -- \u2500\u2500 Step 2: The \"\u22121\" integral vanishes \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\r\n  -- \u222b\u2080^\u03c0 cos(nM) dM = [sin(nM)\/n]\u2080^\u03c0 = 0  (integer n \u2265 1)\r\n  -- Mathlib: integral_cos, Real.sin_nat_mul_pi\r\n  have cos_int_zero :\r\n      \u222b M in (0 : \u211d)..\u03c0, cos (\u2191n * M) = 0 := by\r\n    sorry\r\n\r\n  have step2 :\r\n      \u222b M in (0 : \u211d)..\u03c0, (deriv (eccAnom e) M - 1) * cos (\u2191n * M)\r\n      = \u222b M in (0 : \u211d)..\u03c0, deriv (eccAnom e) M * cos (\u2191n * M) := by\r\n    have key : \u2200 M : \u211d, (deriv (eccAnom e) M - 1) * cos (\u2191n * M)\r\n                       = deriv (eccAnom e) M * cos (\u2191n * M) - cos (\u2191n * M) := by\r\n      intro M; ring\r\n    simp_rw [key]\r\n    rw [intervalIntegral.integral_sub _ _]\r\n    \u00b7 rw [cos_int_zero, sub_zero]\r\n    \u00b7 -- IntervalIntegrable (deriv (eccAnom e) \u00b7 cos(n\u00b7\u00b7))\r\n      sorry\r\n    \u00b7 exact (continuous_cos.comp (continuous_const.mul continuous_id')).intervalIntegrable 0 \u03c0\r\n\r\n  -- \u2500\u2500 Step 3: Change of variable M \u21a6 E via Kepler's equation \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\r\n  -- Under M = E \u2212 e\u00b7sin E:  E'(M) dM = dE,  cos(nM) = cos(nE \u2212 ne\u00b7sinE)\r\n  -- Mathlib: MeasureTheory.integral_image_eq_integral_abs_deriv_smul\r\n  --       or intervalIntegral.integral_comp_deriv\r\n  have step3 :\r\n      \u222b M in (0 : \u211d)..\u03c0, deriv (eccAnom e) M * cos (\u2191n * M)\r\n      = \u222b E in (0 : \u211d)..\u03c0, cos (\u2191n * E - \u2191n * e * sin E) := by\r\n    sorry\r\n\r\n  -- \u2500\u2500 Step 4: Recognise the Bessel integral \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\r\n  -- (1\/\u03c0)\u222b\u2080^\u03c0 cos(nE \u2212 ne\u00b7sinE)dE = J_n(ne)  by definition.\r\n  have step4 :\r\n      (1 \/ \u03c0) * \u222b E in (0 : \u211d)..\u03c0, cos (\u2191n * E - \u2191n * e * sin E)\r\n      = besselJ n (\u2191n * e) := by\r\n    simp only [besselJ]\r\n\r\n  -- \u2500\u2500 Combine all steps \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\r\n  have hpi : \u03c0 > 0 := Real.pi_pos\r\n  have hn_pos : (0 : \u211d) < \u2191n := by exact_mod_cast (show 0 < n by omega)\r\n  rw [step1, step2, step3]\r\n  -- Goal: (2\/\u03c0)\u00b7(1\/n)\u00b7\u222bcos = (2\/n)\u00b7(1\/\u03c0)\u00b7\u222bcos\r\n  -- Both sides equal (2\/(n\u03c0))\u00b7\u222bcos; show by ring after moving (1\/\u03c0) inside.\r\n  have factored : (2 \/ \u03c0) * (1 \/ \u2191n) *\r\n      \u222b E in (0 : \u211d)..\u03c0, cos (\u2191n * E - \u2191n * e * sin E)\r\n      = (2 \/ \u2191n) * ((1 \/ \u03c0) *\r\n      \u222b E in (0 : \u211d)..\u03c0, cos (\u2191n * E - \u2191n * e * sin E)) := by\r\n    ring\r\n  rw [factored, step4]\r\n\r\nend\r\n\r\n\/-! ---------------------------------------------------------------\r\n    \u00a76  Summary of proof obligations\r\n\r\n    Axioms (all standard consequences of the Inverse Function Theorem\r\n    applied to the smooth, strictly increasing map  keplerMap e):\r\n      \u2022 eccAnom           \u2014 existence of the smooth inverse\r\n      \u2022 eccAnom_kepler    \u2014 it satisfies Kepler's equation\r\n      \u2022 eccAnom_hasDerivAt \u2014 C\u00b9 with derivative 1\/(1 \u2212 e\u00b7cos E)\r\n      \u2022 eccAnom_zero      \u2014 E(0) = 0\r\n      \u2022 eccAnom_pi        \u2014 E(\u03c0) = \u03c0\r\n\r\n    `sorry`s (each closes with a standard Mathlib lemma):\r\n      \u2022 step1        IBP via\r\n                      intervalIntegral.integral_deriv_mul_eq_sub_of_hasDerivAt\r\n      \u2022 cos_int_zero  \u222b\u2080^\u03c0 cos(nM)dM = 0 via\r\n                      integral_cos + Real.sin_nat_mul_pi\r\n      \u2022 step2 integrability  IntervalIntegrable for deriv(eccAnom e)\u00b7cos(n\u00b7)\r\n      \u2022 step3        Change of variables via\r\n                      MeasureTheory.integral_image_eq_integral_abs_deriv_smul\r\n\r\n    \u00a77  Finding minimal imports\r\n\r\n    Once the file builds cleanly, add at the bottom:\r\n\r\n        #min_imports\r\n\r\n    and the Lean Infoview will report the exact minimal import list\r\n    for the version of Mathlib you have installed.\r\n    --------------------------------------------------------------- -\/\r\n<\/pre>\n
[1] I used the default setting: Sonnet 4.6, low effort, “thinking” turned off.<\/p>\n","protected":false},"excerpt":{"rendered":"
I ran an experiment today to see whether Claude [1] 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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[202,324],"class_list":["post-247144","post","type-post","status-publish","format-standard","hentry","category-math","tag-artificial-intelligence","tag-formal-methods"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247144","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247144"}],"version-history":[{"count":3,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247144\/revisions"}],"predecessor-version":[{"id":247153,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247144\/revisions\/247153"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247144"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247144"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247144"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247139,"date":"2026-06-10T08:34:31","date_gmt":"2026-06-10T13:34:31","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247139"},"modified":"2026-06-10T13:01:56","modified_gmt":"2026-06-10T18:01:56","slug":"pulling-on-a-thread","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/06\/10\/pulling-on-a-thread\/","title":{"rendered":"Pulling on a thread"},"content":{"rendered":"
Often there’s a thread running through a sequence of my posts. Sometimes I make this explicit and sometimes I don’t.<\/p>\n
The latest thread started with this post<\/a> commenting on a tweet that observed that<\/p>\n
exp(\u2212x<\/em>\u00b2) \u2248 (1 + cos(sin(x<\/em>) + x<\/em>))\/2.<\/p>\n
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<\/a> explaining that it’s not that simple.<\/p>\n
The series for the left hand side alternates and converges very slowly, which lead to the post on naively summing<\/a> an alternating series.<\/p>\n
The series for the right hand side lead to this point on partitions over permutations<\/a>.<\/p>\n
Integrating the right hand side lead to this post<\/a> on how the simplest numerical integration rule works shockingly well on some problems.<\/p>\n
The exact value of the integral turns out to be given by a Bessel function<\/a>, details given in this post<\/a>.<\/p>\n
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<\/a> and on series acceleration<\/a>.<\/p>\n
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.<\/p>\n","protected":false},"excerpt":{"rendered":"
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(\u2212x\u00b2) \u2248 (1 + cos(sin(x) + x))\/2. Some people said online that that the approximation is simply due to the first […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[],"class_list":["post-247139","post","type-post","status-publish","format-standard","hentry","category-math"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247139","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247139"}],"version-history":[{"count":4,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247139\/revisions"}],"predecessor-version":[{"id":247143,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247139\/revisions\/247143"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247139"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247139"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247139"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247135,"date":"2026-06-07T15:14:45","date_gmt":"2026-06-07T20:14:45","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247135"},"modified":"2026-06-07T15:14:45","modified_gmt":"2026-06-07T20:14:45","slug":"aitkin-acceleration-kepler","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/06\/07\/aitkin-acceleration-kepler\/","title":{"rendered":"Aitken acceleration before Aitken"},"content":{"rendered":"
Kepler solved his eponymous equation<\/p>\n
M<\/em> =\u00a0E<\/em> \u2212\u00a0e<\/em> sin(E<\/em>)<\/p>\n
by finding a fixed point of<\/p>\n
E<\/em> = M<\/em> + e<\/em> sin(E<\/em>).<\/p>\n
So guess a value of\u00a0E<\/em> and stick it into the right hand side. Then plug that value into the right hand side again. Kepler said a couple iterations should be enough. And a couple iterations\u00a0are<\/em> enough if the eccentricity e<\/em> is small and you don’t need much accuracy.<\/p>\n
The rate of convergence is determined by\u00a0e<\/em>. Kepler implicitly had in mind small values of\u00a0e<\/em> because he wasn’t aware of anything orbiting the sun in a highly elliptical orbit. Here’s an example with eccentricity 0.05, about the eccentricity of the orbits of Jupiter and Saturn.<\/p>\n
from math import sin\r\n\r\nM, E, e = 1, 1, 0.05\r\nfor _ in range(5):\r\n    E = M + e*sin(E)\r\n    residual = M - (E - e*sin(E))\r\n    print(residual)\r\n<\/pre>\n
The residual after just two iterations is 2.77 \u00d7 10\u22125<\/sup>. If you change e<\/em> to 0.2, the eccentricity of Mercury’s orbit, it takes three iterations to get comparable accuracy. Mercury has the most eccentric orbit of any object Kepler would have known about.<\/p>\n
Now suppose you’d like to solve for E<\/em> when M<\/em> = 1 for Halley’s comet, and you’d like an error of less than 10\u22128<\/sup>. Now you need 16 iterations.<\/p>\n
C. F. W. Peters discovered a faster algorithm in 1891.<\/p>\n
E<\/em>1<\/sub> = M<\/i>\u00a0+ e<\/em> sin(E<\/em>0<\/sub>)
\nE<\/em>2<\/sub> = M<\/i>\u00a0+ e<\/em> sin(E<\/em>1<\/sub>)
\nE<\/em>3<\/sub> = (E<\/em>2<\/sub> E<\/em>0<\/sub> \u2212 E<\/em>1<\/sub>\u00b2)\/(E<\/em>2<\/sub> \u2212 2E<\/em>1<\/sub> + E<\/em>0<\/sub>)<\/p>\n
Let’s look at the results of doing three iterations of Peters’ method for Halley’s comet.<\/p>\n
M, E0, e = 1, 1, 0.967\r\nfor _ in range(3):\r\n    E1 = M + e*sin(E0)\r\n    E2 = M + e*sin(E1)\r\n    E3 = (E2*E0 - E1**2)\/(E2 - 2*E1 + E0)\r\n    residual = M - (E - e*sin(E3))\r\n    print(residual)\r\n    E0 = E3 # for next iteration\r\n<\/pre>\n
This gives a residual of \u22127.23 \u00d7 10\u221210<\/sup>. Each iteration of Peters’ method requires a little more than twice as much work as an iteration of Kepler’s method, but 3 iterations of Peters’ method accomplished more than 16 iterations of Kepler’s method.<\/p>\n
Peters’ algorithm from 1891 was a special case of Alexander Aitken’s series acceleration method published in 1926.<\/p>\n
Related posts<\/h2>\n