Code to generate the pi bound witnesses

Stream: Is there code for X?

[Geoffrey Irving]

I want docs#pi_lt_3141593, but accurate to a bit over 20 decimal digits. However, the file in question seems to lack the code to generate the witnesses to pass to pi_upper_bound and pi_lower_bound. Is that code somewhere else?

[Ruben Van de Velde]

Dates back to https://github.com/leanprover-community/mathlib3/pull/2641. Cc @Mario Carneiro @Floris van Doorn

(Fun fact: in 2021, data.real.pi was the last file to finish compiling in mathlib)

[Daniel Weber]

from sympy import *

a = ceiling(pi * 10**20) / 10**20

v = 0
n = 0
while True:
    ov = sqrt(2+v)
    v = Rational(N(sqrt(2+v) - Integer(10)**-50, 200)).limit_denominator(10**30)
    assert v < ov
    n += 1
    print(v, end=', ')
    if sqrt(2 - v) <= (a - 1 / 4**Integer(n)) / 2 ** (n + 1):
        print('Found', n)
        exit()

Works, and gives

import Mathlib.Data.Real.Pi.Bounds

open Real

theorem pi_lt_314159265358979323847 : π < 3.14159265358979323847 := by
    pi_upper_bound [996006945974988813196889220419/704283265607817876268342353805, 1528727289839560965204682778132/827341247448235339825115450685, 913292375055240730430733739016/465592415232699379941083095141, 1645499627344295663165371448860/826730748193197246562598642091, 1843770496365204984304341838590/922997038538017719302356904659, 1801550067486787619910104613255/901046412076563538906954658887, 1033174046160017800474497982984/516625924061941016458668416829, 1974620922996075539636083349702/987329047688340185766875596479, 1703055998465468187256469972669/851532006704509629847652593643, 1872204633070305500751814237511/936103417906035924020944339290, 417873193942459288085664389516/208936658427256152595763823075, 682652985453241103432668703344/341326517825821018406534464423, 133449833823779055074913079027/66724918138531778002005419822, 1802627778656345399470040951302/901313893470514631525924232247, 1703298290358005799669120926811/851649146157524921098712387678, 1165765145162534233095510819562/582882572748696226264176946563, 171409596335374790028974750983/85704798173841927131273528183, 1916415282813016720770367166539/958207641423710780288743319187, 331175835363368698627932147613/165587917682427537107454982177, 1842175040313531583912899759640/921087520157799292786850855883, 493760775740610311406049235446/246880387870374408374537871617, 1639492360311230529633260363591/819746180155672751777652349511, 1743864522386087812132460629377/871932261193059192729689867963, 24351885390057283142606021981/12175942695028694938275012643, 1393660753917060167524729859317/696830376958530847311606012861, 1199992079080594763280056747597/599996039540297546000874209594, 356242218080878554415612856627/178121109040439289406277007189, 1258407911404392877131542115857/629203955702196449338393661579, 1928835914023314404795982456526/964417957011657206525949495577, 1182695499413497707067870088277/591347749706748854166715368922, 1909465313407151114552936023714/954732656703575557531874429185, 1961572162124030319678167149025/980786081062015159904677606876, 931266043180727152209784389970/465633021590363576112677467451, 1670978975783135072358731367231/835489487891567536182857979589]

For lower bounds you can use

from sympy import *

a = floor(pi * 10**20) / 10**20

v = 0
n = 0
while True:
    ov = sqrt(2+v)
    v = Rational(N(sqrt(2+v) + Integer(10)**-50, 200)).limit_denominator(10**30)
    assert v > ov
    n += 1
    print(v, end=', ')
    if sqrt(2 - v) >= (a) / 2 ** (n + 1):
        print('Found', n)
        exit()

[Geoffrey Irving]

Thank you!

[Geoffrey Irving]

Hmm, is there a typo in the docstring for Real.«tacticPi_upper_bound[_,,]»? I feel like that last ≥ should be ≤.

Create a proof of π < a for a fixed rational number a, given a witness, which is a sequence of rational numbers √2 < r 1 < r 2 < ... < r n < 2 satisfying the property that √(2 + r i) ≥ r(i+1), where r 0 = 0 and √(2 - r n) ≥ (a - 1/4^n) / 2^(n+1).

[Yury G. Kudryashov]

Would it be better to use a series that converges to π faster? E.g., there are various expressions with atans, some of them converge pretty fast.

[Mario Carneiro]

Geoffrey Irving said:

I want docs#pi_lt_3141593, but accurate to a bit over 20 decimal digits. However, the file in question seems to lack the code to generate the witnesses to pass to pi_upper_bound and pi_lower_bound. Is that code somewhere else?

I wrote some adhoc code in mathematica to optimize the witnesses back in the day. I doubt I still have it but I can try to reconstruct it

[Geoffrey Irving]

Daniel’s code above works fine for my purposes (and the resulting proof is still fast). So quite low priority!

[Ruben Van de Velde]

If you pr that bound, can you add the script in a comment?

[Eric Wieser]

What's the obstruction in translating @Daniel Weber's example to lean; is it just the limit_denominator?

[Geoffrey Irving]

N and sqrt are obstacles as well.

[Geoffrey Irving]

And pi, of course: one would need a separate (unverified) series for pi, though that’s easy.

[Daniel Weber]

For Lean it would probably be better to replace limit_denominator by Stern–Brocot tree search, and then you only need to compare to the square root

[Mario Carneiro]

Here's an optimization of the proof in this case, containing a few new tricks:

theorem pi_lt_3141593' : π < 3.14159265358979323847 := by
  pi_upper_bound [
    215157040700/152139002499, 936715022285/506946517009, 1760670193473/897581880893,
    2-6049918861/628200981455, 2-8543385003/3546315642356, 2-2687504973/4461606579043,
    2-1443277808/9583752057175, 2-546886849/14525765179168, 2-650597193/69121426717657,
    2-199969519/84981432264454, 2-226282901/384655467333100, 2-60729699/412934601558121,
    2-25101251/682708800188252, 2-7156464/778571703825145, 2-7524725/3274543383827551,
    2-4663362/8117442793616861, 2-1913009/13319781840326041, 2-115805/3225279830894912,
    2-708749/78957345705688293, 2-131255/58489233342660393, 2-101921/181670219085488669,
    2-44784/319302953916238627, 2-82141/2342610212364552264, 2-4609/525783249231842696,
    2-4567/2083967975041722089, 2-2273/4148770928197796067, 2-563/4110440884426500846,
    2-784/22895812812720260289, 2-1717/200571992854289218531, 2-368/171952226838388893139,
    2-149/278487845640434185590, 2-207/1547570041545500037992, 2-20/598094702046570062987,
    2-7/837332582865198088180]

[Mario Carneiro]

Here's the Mathematica code that produced it, after some cleanup:

calc[a_, Iters -> n_, Rounding -> extra_, Precision -> prec_] :=
  Module[{r0, r, r2, diff},
   On[Assert];
   r0 = 2 - ((a - 1/4^n)/2^(n + 1))^2;
   r = Log[2 - NestList[((#)^2 - 2) &, N[r0, prec], n - 1]];
   If[r[[-1]] < Log[2 - Sqrt[2]], Return["insufficient iterations"]];
   diff = (r[[-1]] - Log[2 - Sqrt[2]])/(Length[r] + 1);
   r2 = Log[
     Rationalize[Exp[#], extra (Exp[# + diff] - Exp[#])] &
        /@ (r - diff Range[1, Length[r]])];
   Assert[2 - Exp@r2[[1]] >= r0];
   Assert[
    And @@ Table[
      Sqrt@(4 - Exp@r2[[i + 1]]) >= (2 - Exp@r2[[i]]), {i, 1,
       Length[r2] - 1}]];
   Assert[Exp@r2[[-1]] >= 2 - Sqrt[2]];
   With[{s1 = ToString@InputForm[2 - #], s2 = ToString@InputForm[#]},
      If[StringLength[s1] <= StringLength[s2] + 2, s1,
       "2-" <> s2]] & /@ Exp@Reverse@r2];
calc[314159265358979323847/100000000000000000000, Iters -> 34,
 Rounding -> .5, Precision -> 46]

[Mario Carneiro]

I also had to slightly modify pi_upper_bound to allow rational expressions:

/-- Create a proof of `π < a` for a fixed rational number `a`, given a witness, which is a
sequence of rational numbers `√2 < r 1 < r 2 < ... < r n < 2` satisfying the property that
`√(2 + r i) ≥ r(i+1)`, where `r 0 = 0` and `√(2 - r n) ≥ (a - 1/4^n) / 2^(n+1)`. -/
elab "pi_upper_bound " "[" l:term,* "]" : tactic => do
  have els := l.getElems
  let n := quote els.size
  evalTactic (← `(tactic| apply pi_upper_bound_start $n))
  for l in els do
    let Q := .const ``_root_.Rat []
    let {num, den, ..} ← unsafe Meta.evalExpr ℚ Q (← Term.elabTermAndSynthesize l (some Q))
    evalTactic (← `(tactic| apply sqrtTwoAddSeries_step_down $(quote num.toNat) $(quote den)))
  evalTactic (← `(tactic| simp [sqrtTwoAddSeries]))
  allGoals do
    evalTactic (← `(tactic| (first | norm_num1 | done)))

[Eric Wieser]

Could docs#Real.sqrtTwoAddSeries_step_down just take the Rat directly?

[Mario Carneiro]

the Rat processing is done at compile time here, the resulting proof uses only nats

[Mario Carneiro]

and the reason for this is so that it can get cross multiplied away in the subgoal which is proven by kernel computation

[Eric Wieser]

Where's that kernel computation happening in your elab above?

[Mario Carneiro]

well, this is not a very high performance implementation, it got dumbed down a bit during the port. But it's still doing kernel computation via norm_num1

[Eric Wieser]

Oh, I thought that norm_num avoided kernel computation and just dutifully chained a lot of lemmas together

[Yury G. Kudryashov]

norm_num uses rfl to prove various steps.

[Yury G. Kudryashov]

It tries to replace some heavy rfls with something easier to check, but there are kernel computations involved anyway.

[Mario Carneiro]

in lean 3 it was just lemma application, but lean 4 norm_num knows about kernel accelerated nat ops and uses them whenever possible

[Mario Carneiro]

#16930 adds all of the above to the pi bounds file. @Geoffrey Irving can you clarify if the bounds are good enough for your application, since you said "a bit over 20 decimal digits"?

[Geoffrey Irving]

Confirmed, those bounds are good enough for me. When I said "a bit over 20 decimal digits", I had forgotten that it is not the case that π < 1 (I needed relative error).

[Geoffrey Irving]

And you don't even need 20 digits to know that!