happy pi day!

Stream: general

[Kenny Lau]

happy pi day!

[Kevin Buzzard]

https://www.bbc.co.uk/news/technology-47524760

[Kevin Buzzard]

But did they formally verify it? And if not, should we believe the calculations? They should be reproduced using a completely different system.

[Moses Schönfinkel]

They use https://en.wikipedia.org/wiki/Bellard%27s_formula to verify digits they compute with https://en.wikipedia.org/wiki/Chudnovsky_algorithm

[Moses Schönfinkel]

Also it would appear that they are simply using Alexander Yee's pi cruncher to simply run their computation, this is its homepage: http://www.numberworld.org/y-cruncher/

[Koundinya Vajjha]

But did they formally verify it? And if not, should we believe the calculations? They should be reproduced using a completely different system.

https://arxiv.org/abs/1709.01743

[Moses Schönfinkel]

There's also an extra verification algorithm they employ it seems, so you can trust it about as much as you can trust the (un-verified C++) implementation of Lean that runs on un-verified hardware.

[Kenny Lau]

Can we prove that pi starts with 3.14? i.e. 3.14 <= pi < 3.15?

[Kenny Lau]

I think a Chinese mathematician 1000 years ago used 3•2^n-gons to estimate pi up to 7 digits?

[Chris Hughes]

The easiest way using current mathlib, is to just use some lemma I proved about approximating e^x for x le 1, as the first n terms of the power series. I think tactics to do this stuff would be good.

[Kevin Buzzard]

Can we prove that pi starts with 3.14? i.e. 3.14 <= pi < 3.15?

Do we have any theorem of the form "pi is the following infinite sum" or "pi is the limit of the following sequence"?

[Kevin Buzzard]

Maybe today is a good day to talk about pi. We know cos(pi/2)=0 so we can probably deduce any reasonable trig assertion about pi. Now the derivative of arctan is 1/(1+x^2) and...oh dear now I want to integrate term by term to deduce the power series for arctan, and then evaluate at 1 to deduce pi/4=1-1/3+1/5-1/7+..., for which one can get some reasonable approximations to pi. But we still need derivative of a power series :-/

[Kenny Lau]

that power series is notoriously slow to converge

[Kevin Buzzard]

I'm well aware of that. Conversely it's an alternating series so it provides us with a rigorous proof that, for example, pi>1-1/3+1/5-1/7+1/9-1/11 and pi<1-1/3+1/5-1/7+1/9-1/11+1/13.

[Kevin Buzzard]

Which is more than we know currently in mathlib.

[Kevin Buzzard]

The 2^n-gon trick converges much more quickly and we have the trig in place, but it relies on the fact that pi = circumference over diameter, which we can't formalise yet.

[Kevin Buzzard]

Claim: pi > 3. Proof: draw a regular hexagon and chop it into 6 triangles. Then draw a circle round the hexagon.

[Kevin Buzzard]

Circumference of circle > circumference of hexagon = 6, diameter = 2.

[Kevin Buzzard]

Can we turn this idea into a mathlib proof that pi > 3? If so then the polygon trick starts to kick in.

[Abhimanyu Pallavi Sudhir]

You'll need pi as an integral, correct?

[Kevin Buzzard]

I don't know. We've defined pi "in terms of trig", and we're "doing trig" in some vague way above -- but it looks to me like we're somehow doing a different trig.

[Kenny Lau]

just define it as the limit of the 3*2^n-gons

[Abhimanyu Pallavi Sudhir]

I think the easiest way to prove the circumference of that limit is pi still uses Riemann sums, though.

[Kevin Buzzard]

ha ha

[Kevin Buzzard]

there's a problem with the polygon approach

[Kevin Buzzard]

you end up defining 2pi as the limit, as n tends to infinity, of $n\sin(2\pi/n)$

[Kevin Buzzard]

(that's using n-gons; using 2^n-gons just change n to 2^n)

[Kevin Buzzard]

OTOH if one could prove that $n\sin(2\pi/n)$ was increasing

[Kevin Buzzard]

and that the derivative of sin was cos

[Kevin Buzzard]

then setting n=6 you get pi>3

[Reid Barton]

For lower bounds, I guess you just need to prove $\sin x < x$ for $x > 0$ which seems feasible

[Kevin Buzzard]

aah so this is the point; we can compute sin(2 pi / (3*2^m)) as a solution to a quadratic

[Kevin Buzzard]

For lower bounds, I guess you just need to prove $\sin x < x$ for $x > 0$ which seems feasible

yeah, just draw the graph :P

[Kevin Buzzard]

that crazy quote bug again

[Kevin Buzzard]

sin(x)<x should be easy by that alternating series thing

[Kevin Buzzard]

Maybe Chris even proved sin(x)>x-x^3/3 which might even have been used to prove the current known bounds for pi, which are something like 2<pi<4

[Kevin Buzzard]

good to see that google are on 33 trillion digits and we're still working on the 0th one :-)

[Kevin Buzzard]

the thing is, when we get it, we'll have a proof

[Kevin Buzzard]

not just some blather about cloud computing

[Kevin Buzzard]

It's Xena tonight, maybe we should have a pi day celebration and try to prove pi > 3

[Kenny Lau]

firstly how do you define pi

[Kevin Buzzard]

it's in mathlib

[Kevin Buzzard]

so it is defined by official decree already

[Kevin Buzzard]

I think they went for 22/7 but I can't remember

[Kenny Lau]

it's in mathlib?

[Kevin Buzzard]

no wait, that was something else.

[Kevin Buzzard]

Isn't it in mathlib??

[Kevin Buzzard]

yes it's there

[Kevin Buzzard]

It's in analysis.exponential

[Kevin Buzzard]

local notation π := 22 / 7

[Kevin Buzzard]

no wait

[Kenny Lau]

oh right they proved sin(pi) = 0 right

[Kevin Buzzard]

local notation π := real.pi

[Kenny Lau]

oh

[Kenny Lau]

what's the definition?

[Abhimanyu Pallavi Sudhir]

Two times the zero of cos x between 0 and 2. I don't see why they didn't just prove exp periodic and define it as the half-period. Maybe that's hard without knowing what pi is, but it feels more natural.

[Kevin Buzzard]

I think they were pressed for time and wanted to get it out for my birthday :-)

[Kevin Buzzard]

I gave a talk yesterday to 100 people at my university interested in teaching, and told them the story of how I got pi for my birthday, and it got the biggest laugh of the talk.

[Kevin Buzzard]

You weren't around at that point Abhi, so you couldn't keep them in order :-)

[Kevin Buzzard]

noncomputable def pi : ℝ := 2 * classical.some exists_cos_eq_zero

[Kevin Buzzard]

holy crap, it could be twice any zero at all!

[Kevin Buzzard]

Maybe it's not less than 4 after all!

[Kevin Buzzard]

pesky axiom of choice

[Kevin Buzzard]

lemma exists_cos_eq_zero : ∃ x, 1 ≤ x ∧ x ≤ 2 ∧ cos x = 0 :=

[Kevin Buzzard]

Proof uses cos_one_pos and cos_two_neg and IVT

[Kevin Buzzard]

so indeed currently we should be able to prove that pi is between 2 and 4

[Kevin Buzzard]

but for proving it's bigger than 3 it does not suffice to prove that cos(1.5) is positive, because there might be some other zero of cos in the range (1,1.5).

[Kevin Buzzard]

It would be good to know that cos was decreasing. That would follow from MVT perhaps, if we knew the derivative of cos was sin, but we might be a long way from that.

[Kevin Buzzard]

Do we know the derivative of exp is exp yet?

[Abhimanyu Pallavi Sudhir]

I don't. Maybe Jeremy Avigad does (if he's specialised his derivative for basic functions).

[Reid Barton]

I think we do know cos is decreasing, or at least I saw something similar when looking through the file for sin x < x

[Reid Barton]

cos_lt_cos_of_nonneg_of_le_pi_div_two

[Reid Barton]

or actually cos_lt_cos_of_nonneg_of_le_pi

[Reid Barton]

It's proved using angle addition formulas

[Floris van Doorn]

If we know that sin x < x here is a proof sketch to show that pi > 3:

The doubling formula cos(2*x) = 2*cos^2(x) - 1 has been formalized, so we can compute cos^2(pi/4) and cos^2(pi/8). Now we can also compute sin^2(pi/8). This turns out to be 1/2 - 1/sqrt(8). Now we know that pi > 8sin(pi/8) = 8sqrt(1/2-1/sqrt(8)) > 3.

[Kenny Lau]

nice!

[Chris Hughes]

I can imagine this being a long annoying proof.

[Reid Barton]

sin_bound should be enough to prove sin x < x as well

[Floris van Doorn]

I'm starting now to compute sin(pi/8).

[Kevin Buzzard]

Happy pi day :-)

[Kevin Buzzard]

Chris this proof doesn't look too bad to me at all. I was manipulating 1+sqrt(5)/2 without too much trouble last week when we were doing some analysis example sheet question.

[Kevin Buzzard]

Or is there something else you're worried about?

[Floris van Doorn]

How do I prove (2 : ℂ) ≠ 0?

[Kevin Buzzard]

These things should go in pi.lean

[Kevin Buzzard]

norm_num proves 2 ne 0

[Kevin Buzzard]

or if it doesn't then take real parts

[Kevin Buzzard]

because it will definitely prove (2 : real) \ne 0

[Chris Hughes]

two_ne_zero'

[Kevin Buzzard]

rofl

[Chris Hughes]

It needs the prime, because I think the other version is for linear ordered things.

[Floris van Doorn]

yes, I found the other one. Thanks

[Kevin Buzzard]

linear ordered -> char 0

[Chris Hughes]

The linear order version is in core I think.

[Kevin Buzzard]

https://leanprover.zulipchat.com/#narrow/stream/116395-maths/topic/golden.20ratio.20calculation

[Kevin Buzzard]

@Floris van Doorn there's an example of basic manipulations with sqrt(5)

[Floris van Doorn]

import analysis.exponential

namespace real
variable (x : ℝ)

lemma sqrt_le_left {x y : ℝ} (hy : 0 ≤ y) : sqrt x ≤ y ↔ x ≤ y ^ 2 :=
begin
  rw [mul_self_le_mul_self_iff (sqrt_nonneg _) hy, pow_two],
  cases le_total 0 x with hx hx,
  { rw [mul_self_sqrt hx] },
  { have h1 : 0 ≤ y * y := mul_nonneg hy hy,
    have h2 : x ≤ y * y := le_trans hx h1,
    simp [sqrt_eq_zero_of_nonpos, hx, h1, h2] }
end

lemma le_sqrt {x y : ℝ} (hx : 0 ≤ x) (hy : 0 ≤ y) : x ≤ sqrt y ↔ x ^ 2 ≤ y :=
by rw [mul_self_le_mul_self_iff hx (sqrt_nonneg _), pow_two, mul_self_sqrt hy]

lemma div_lt_self {x y : ℝ} (hx : x > 0) (hy : y > 1) : x / y < x :=
by simpa using div_lt_div' (le_refl x) hy hx zero_lt_one

lemma cos_square : cos x ^ 2 = (1 + cos (2 * x)) / 2 :=
by simp [cos_two_mul, mul_div_cancel_left, two_ne_zero]

lemma sin_square : sin x ^ 2 = 1 - cos x ^ 2 :=
by { rw [←sin_pow_two_add_cos_pow_two x], simp }

lemma cos_square_pi_div_four : cos (pi / 4) ^ 2 = 1 / 2 :=
by { have : 2 * (pi / 4) = pi / 2, ring, simp [cos_square, this] }

lemma cos_pi_div_four : cos (pi / 4) = sqrt (1 / 2) :=
begin
  symmetry, rw [sqrt_eq_iff_sqr_eq, cos_square_pi_div_four], norm_num,
  apply le_of_lt, apply cos_pos_of_neg_pi_div_two_lt_of_lt_pi_div_two,
  { transitivity (0 : ℝ), rw neg_lt_zero, apply pi_div_two_pos,
    exact div_pos pi_pos (by norm_num) },
  apply div_lt_div' (le_refl pi) _ pi_pos _,
  norm_num, norm_num
end

/- why is division_ring.inv_inj in the file "field"?-/
lemma cos_square_pi_div_eight : cos (pi / 8) ^ 2 = 1 / 2 + sqrt (1 / 8) :=
begin
  have h1 : 2 * (pi / 8) = pi / 4, ring,
  have h2 : sqrt 4 = 2, rw [sqrt_eq_iff_sqr_eq]; norm_num,
  have h3 : 2 * sqrt 2 ≠ 0, norm_num,
  have h4 : sqrt 8 ≠ 0, norm_num,
  simp [cos_square, h1, cos_pi_div_four],
  rw [←div_add_div_same], simp [div_eq_mul_inv, (mul_inv' _ _).symm, division_ring.inv_inj, h3, h4],
  transitivity sqrt 4 * sqrt 2, rw [h2],
  rw [←sqrt_mul]; norm_num
end

lemma cos_pi_div_eight : cos (pi / 8) = sqrt (1 / 2 + sqrt (1 / 8)) :=
begin
  symmetry, rw [sqrt_eq_iff_sqr_eq, cos_square_pi_div_eight],
  apply add_nonneg, norm_num, apply sqrt_nonneg,
  apply le_of_lt, apply cos_pos_of_neg_pi_div_two_lt_of_lt_pi_div_two,
  { transitivity (0 : ℝ), rw neg_lt_zero, apply pi_div_two_pos,
    exact div_pos pi_pos (by norm_num) },
  apply div_lt_div' (le_refl pi) _ pi_pos _,
  norm_num, norm_num
end

lemma sin_square_pi_div_eight : sin (pi / 8) ^ 2 = 1 / 2 - sqrt (1 / 8) :=
by { simp [sin_square, cos_square_pi_div_eight], rw [←add_assoc], congr' 1, norm_num }

lemma sin_pi_div_eight : sin (pi / 8) = sqrt (1 / 2 - sqrt (1 / 8)) :=
begin
  symmetry, rw [sqrt_eq_iff_sqr_eq, sin_square_pi_div_eight],
  { rw [sub_nonneg, sqrt_le_left], norm_num, norm_num },
  apply le_of_lt, apply sin_pos_of_pos_of_lt_pi,
  { exact div_pos pi_pos (by norm_num) },
  exact div_lt_self pi_pos (by norm_num)
end

lemma sin_lt {x : ℝ} (h : 0 < x) : sin x < x :=
begin
  cases le_or_gt x 1 with h' h',
  { have hx : abs x = x := abs_of_nonneg (le_of_lt h),
    have : abs x ≤ 1, rwa [hx],
    have := sin_bound this, rw [abs_le] at this,
    have := this.2, rw [sub_le_iff_le_add', hx] at this,
    apply lt_of_le_of_lt this, rw [sub_add], apply lt_of_lt_of_le _ (le_of_eq (sub_zero x)),
    apply sub_lt_sub_left, rw sub_pos, apply mul_lt_mul',
    { rw [pow_succ x 3], refine le_trans _ (le_of_eq (one_mul _)),
      rw mul_le_mul_right, exact h', apply pow_pos h },
    norm_num, norm_num, apply pow_pos h },
  exact lt_of_le_of_lt (sin_le_one x) h'
end

lemma pt_gt_three : pi > 3 :=
begin
  have : pi > sin (pi / 8) * 8,
  { apply mul_lt_of_lt_div, norm_num, apply sin_lt, apply div_pos pi_pos, norm_num },
  apply lt_of_le_of_lt _ this, rw [sin_pi_div_eight],
  apply le_mul_of_div_le, norm_num,
  rw [le_sqrt, le_sub_iff_add_le, ←le_sub_iff_add_le', sqrt_le_left],
  norm_num, norm_num, norm_num,
  rw [sub_nonneg, sqrt_le_left], norm_num, norm_num
end

end real

[Kevin Buzzard]

This doesn't compile for me

[Kevin Buzzard]

solve1 tactic failed, focused goal has not been solved
state:
⊢ pi / 8 ≤ 1

[Kevin Buzzard]

We need pi <= 8

[Kevin Buzzard]

lemma pt_gt_three : pi > 3 :=
begin
  have : pi > sin (pi / 8) * 8,
  { apply mul_lt_of_lt_div, norm_num, apply sin_lt, apply div_pos pi_pos, norm_num,
  -- ⊢ pi / 8 ≤ 1
    sorry
  },
  apply lt_of_le_of_lt _ this, rw [sin_pi_div_eight],
  apply le_mul_of_div_le, norm_num,
  rw [le_sqrt, le_sub_iff_add_le, ←le_sub_iff_add_le', sqrt_le_left],
  norm_num, norm_num, norm_num,
  rw [sub_nonneg, sqrt_le_left], norm_num, norm_num
end

[Floris van Doorn]

oops, I missed that. One minute

[Floris van Doorn]

I updated the previous post, it should compile now.

[Reid Barton]

You can also prove sin x < x for all x > 0 by handling the case x > 1 by sin x <= 1 < x

[Floris van Doorn]

@Reid Barton That is indeed the nicer solution. I was planning to do that, but I was a bit too eager to post it here and forgot :)

[Floris van Doorn]

I updated the code once more.

[Kevin Buzzard]

This should all go in pi.lean, right?

[Kevin Buzzard]

So now we can prove that floor(pi) = 3

[Kevin Buzzard]

So the next question is whether we can prove digit 10 pi 0 = 1

[Kevin Buzzard]

ie floor(pi*10)%10=1

[Johan Commelin]

This should all go in pi.lean, right?

Isn't that about pi types :rolling_on_the_floor_laughing: :cake:

[Floris van Doorn]

Computing sin(pi/16) we can conclude that pi > 3.1, computing sin(pi/64) we can conclude that pi>3.14. Using the other side of sin_bound we can probably also get upper bounds.

[Kevin Buzzard]

Can you get to 32 trillion by the end of pi day?

[Floris van Doorn]

probably not

[Floris van Doorn]

Is there a better way than norm_num to prove the following?

(2 : ℝ) ≤ ((((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2

(the RHS is approx 2.00033). norm_num is struggling on this one, and gets a deterministic timeout.

[Floris van Doorn]

If I can prove inequalities like that, I can approximate pi for a few decimal places, at least from below (I need to still work on approximation from above). The approximation I have of pi is of the form

π ≈ 2 ^ 7 * sqrt (2 - sqrt (2 + sqrt (2 + sqrt (2 + sqrt (2 + sqrt (2 + sqrt 2))))))

This seems to converge at an okay rate: the number of correct digits seems to be linear in the number of square roots.

In particular, I formalized

lemma pt_gt_314 : (2 : ℝ) ≤ ((((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 → pi > 3.14

[Kevin Buzzard]

I can reduce your time-out goal to ⊢ 0 < 3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000

[Kevin Buzzard]

It's not even close

[Kevin Buzzard]

but I'm struggling with the last bit.

[Kevin Buzzard]

oh -- these are reals

[Kevin Buzzard]

import data.real.basic

example : (2 : ℝ) ≤ ((((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 :=
begin
  rw (show 4 + 1 = 5, from rfl),
  rw (show (2 : ℝ) ^ 5 = 32, by norm_num),
  rw (show (157 : ℝ) / 50 / 32 = 157 / 1600, by norm_num),
  rw (show ((157 : ℝ) / 1600) ^ 2 = 24649 / 2560000, by norm_num),
  rw (show (2 : ℝ) - 24649 / 2560000 = 5095351/2560000, by norm_num),
  rw (show ((5095351 : ℝ) /2560000) ^ 2 = 25962601813201/6553600000000, by norm_num),
  -- let's try skipping steps
  rw (show ((((25962601813201 : ℝ) / 6553600000000 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 =
    6806775565993596917798714352237438749189222725565123913061124308255143227446400872965401873904861225601/3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000,
    by norm_num),
  -- worked!
  -- ⊢ 2 ≤
  --  6806775565993596917798714352237438749189222725565123913061124308255143227446400872965401873904861225601 /
  --    3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
  rw le_div_iff,
    norm_num,
  -- ⊢ 0 < 3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
  sorry
end

[Chris Hughes]

norm_num just works for that last goal for me.

example : (2 : ℝ) ≤ ((((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 :=
begin
  rw (show 4 + 1 = 5, from rfl),
  rw (show (2 : ℝ) ^ 5 = 32, by norm_num),
  rw (show (157 : ℝ) / 50 / 32 = 157 / 1600, by norm_num),
  rw (show ((157 : ℝ) / 1600) ^ 2 = 24649 / 2560000, by norm_num),
  rw (show (2 : ℝ) - 24649 / 2560000 = 5095351/2560000, by norm_num),
  rw (show ((5095351 : ℝ) /2560000) ^ 2 = 25962601813201/6553600000000, by norm_num),
  -- let's try skipping steps
  rw (show ((((25962601813201 : ℝ) / 6553600000000 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 =
    6806775565993596917798714352237438749189222725565123913061124308255143227446400872965401873904861225601/3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000,
    by norm_num),
  -- worked!
  -- ⊢ 2 ≤
  --  6806775565993596917798714352237438749189222725565123913061124308255143227446400872965401873904861225601 /
  --    3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
  rw le_div_iff,
  { norm_num },
  { norm_num }
  -- ⊢ 0 < 3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000

end

[Kevin Buzzard]

import data.real.basic

--set_option pp.all true
example : (2 : ℝ) ≤ ((((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 :=
begin
  rw (show 4 + 1 = 5, from rfl),
  rw (show (2 : ℝ) ^ 5 = 32, by norm_num),
  rw (show (157 : ℝ) / 50 / 32 = 157 / 1600, by norm_num),
  rw (show ((157 : ℝ) / 1600) ^ 2 = 24649 / 2560000, by norm_num),
  rw (show (2 : ℝ) - 24649 / 2560000 = 5095351/2560000, by norm_num),
  rw (show ((5095351 : ℝ) /2560000) ^ 2 = 25962601813201/6553600000000, by norm_num),
  -- let's try skipping steps
  rw (show ((((25962601813201 : ℝ) / 6553600000000 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 =
    6806775565993596917798714352237438749189222725565123913061124308255143227446400872965401873904861225601/3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000,
    by norm_num),
  -- worked!
  -- ⊢ 2 ≤
  --  6806775565993596917798714352237438749189222725565123913061124308255143227446400872965401873904861225601 /
  --    3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
  rw le_div_iff,
    norm_num,
  -- ⊢ 0 < 3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
  rw (show (3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000 : ℝ)
           = 40 ^ 64, by norm_num),
  apply pow_pos,
  norm_num,
end

This worked for me.

[Kevin Buzzard]

Oh, I think that I upset Lean with some overflow and then failed to restart, and had problems.

[Kevin Buzzard]

Only 31999999999998 digits to go and we have the record!

[Kenny Lau]

oh boy

[Floris van Doorn]

I have proven

lemma pi_gt_314 : pi > 3.14
lemma pi_lt_315 : pi < 3.15

assuming the following 12 inequalities:

⊢ 2 ≤ ((((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2
⊢ 0 ≤ (((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2
⊢ 0 ≤ ((2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2
⊢ 0 ≤ (2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2
⊢ 0 ≤ 2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2
⊢ 0 < 157 / 50 / 2 ^ (4 + 1)
⊢ ((((2 - ((63 / 20 - 1 / 4 ^ 4) / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2) ^ 2 ≤ 2
⊢ 0 < (((2 - ((63 / 20 - 1 / 4 ^ 4) / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2) ^ 2 - 2
⊢ 0 < ((2 - ((63 / 20 - 1 / 4 ^ 4) / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2) ^ 2 - 2
⊢ 0 < (2 - ((63 / 20 - 1 / 4 ^ 4) / 2 ^ (4 + 1)) ^ 2) ^ 2 - 2
⊢ 0 < 2 - ((63 / 20 - 1 / 4 ^ 4) / 2 ^ (4 + 1)) ^ 2
⊢ 0 ≤ (63 / 20 - 1 / 4 ^ 4) / 2 ^ (4 + 1)

[Kenny Lau]

I've verified all of your inequalities on WolframAlpha

[Mario Carneiro]

Oh boy I've missed a lot

[Mario Carneiro]

For these kinds of numeric calculations you can simplify a lot by decreasing the precision in the middle of the calculation

[Mario Carneiro]

Are all the inequalities implied by one of them?

[Kenny Lau]

the first one is exceptionally tight

[Kenny Lau]

the 2nd to 5th looks like a more general theorem?

[Kenny Lau]

wait are those approximations to sqrt(2)?

[Kenny Lau]

@Floris van Doorn I think it would be better if you tell us more about what's actually going on

[Kenny Lau]

if we want approximations to sqrt(2) maybe we should turn to Pell theory

[Mario Carneiro]

I think they are approximations to 2

[Mario Carneiro]

fixed points of x |-> x^2-2

[Kenny Lau]

interesting...?

[Mario Carneiro]

The first one isn't that tight: it is saying that

   1128227574827648531222203602074520069222725565123913061124308255143227446400872965401873904861225601 /
3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
>=0

[Mario Carneiro]

Notice that there are like 60 digits of precision there and only 4 of them matter

[Mario Carneiro]

so you can vastly simplify the calculation by precision reduction in the middle

[Kenny Lau]

aha

[Mario Carneiro]

π ≈ 2 ^ 7 * sqrt (2 - sqrt (2 + sqrt (2 + sqrt (2 + sqrt (2 + sqrt (2 + sqrt 2))))))

Where did this come from?

[Floris van Doorn]

The file I'm working on is here: https://github.com/fpvandoorn/mathlib/blob/pi/src/data/real/pi.lean

[Kenny Lau]

/- the series sqrt(2+sqrt(2+ ... ))-/
noncomputable def sqrt_two_add_series : ℕ → ℝ
| 0     := 0
| (n+1) := sqrt (2 + sqrt_two_add_series n)

lemma pi_gt_sqrt_two_add_series (n : ℕ) : pi > 2 ^ (n+1) * sqrt (2 - sqrt_two_add_series n) :=

lemma pi_lt_sqrt_two_add_series (n : ℕ) : pi < 2 ^ (n+1) * sqrt (2 - sqrt_two_add_series n) + 1 / 4 ^ n :=

lemma pi_gt_314 : pi > 3.14 :=
begin
  refine lt_of_le_of_lt _ (pi_gt_sqrt_two_add_series 4),

lemma pi_lt_315 : pi < 3.15 :=
begin
  refine lt_of_lt_of_le (pi_lt_sqrt_two_add_series 4) _,

[Kenny Lau]

extract ^

[Kenny Lau]

anyway @Floris van Doorn remarkable work!

[Mario Carneiro]

The initial starting point sqrt_two_add_series 0 is arbitrary? What constraints does it require?

[Floris van Doorn]

The first and 7th inequalities are relatively tight (it compares 2 with 2.0003 and 1.994). The other inequalities are not tight at all, the numbers on the right are all between 1 and 2.
What's going on is that I proved that

2 ^ 5 * sqrt (2 - sqrt (2 + sqrt (2 + sqrt (2 + sqrt (2)))) < π < 2 ^ 5 * sqrt (2 - sqrt (2 + sqrt (2 + sqrt (2 + sqrt (2)))) + 1 / 4 ^ 4

Now we can check the bounds by replacing π by 3.14 and 3.15 and then moving all the square roots of the inequality to the other side. To do that, we need to check that stuff is positive (that gives all the non-tight inequalities), and the final inequality are the two tight ones.

[Mario Carneiro]

If we let sqrt_two_add_series a n be the series starting at a, then my hope is to get an inequality chain based on sqrt_two_add_series a (n+1) <= sqrt_two_add_series b n when a and b are suitably related rational numbers

[Mario Carneiro]

is sqrt_two_add_series a n monotonic in a? I think it is

[Floris van Doorn]

The defining equation of sqrt_two_add_series is

cos (pi / 2 ^ (n+1)) = sqrt_two_add_series n / 2

[Mario Carneiro]

aha, you inverted sin^2+cos^2 = 1

[Floris van Doorn]

yes, the main equality is lemma cos_square : cos x ^ 2 = 1 / 2 + cos (2 * x) / 2

[Kenny Lau]

how do we formalize "precision reduction"? @Mario Carneiro

[Mario Carneiro]

blasty lean users aren't going to like the answer

[Mario Carneiro]

you give intermediate small rational numbers as certificates

[Floris van Doorn]

yes, it might be nice to let sqrt_two_add_series start at an arbitrary number, and yes, sqrt_two_add_series a n is monotone in both a and n.

[Kenny Lau]

fortunately I'm not a blasty lean user

[Mario Carneiro]

Is it true that sqrt_two_add_series a (n+1) = sqrt_two_add_series (sqrt a + 2) n or something like that?

[Mario Carneiro]

sorry playing catchup

[Floris van Doorn]

almost: sqrt_two_add_series a (n+1) = sqrt_two_add_series (sqrt (2 + a)) n

[Mario Carneiro]

Okay, so with monotonicity we get the desired result: if 2+a <= b^2 and b >= 0 then sqrt_two_add_series a (n+1) <= sqrt_two_add_series b n

[Mario Carneiro]

so now you apply this lemma 5 times, and supply 5 small rational numbers that satisfy the inequalities (check in mathematica or something)

[Mario Carneiro]

and those numbers are the certificate

[Floris van Doorn]

that sounds good

[Floris van Doorn]

(of course I hoped for a lazier solution :) )

[Mario Carneiro]

Oh the lazy solution is also possible, you can just tell lean "round to n digits at each stage"

[Mario Carneiro]

and pick n so it works

[Mario Carneiro]

but the really really important part is don't work out all the algebra with the original number, that causes exponential blowup

[Mario Carneiro]

I really like certificate proofs, they show the added power of ITP

[Kenny Lau]

didn't you write a floating point library or something?

[Mario Carneiro]

you can solve NP problems in polynomial time, people! That's big news

[Kenny Lau]

would that be useful?

[Mario Carneiro]

this is why metamath is fast and every tactic based prover is 5 orders of magnitude slower, because certificates matter

[Kenny Lau]

so we still use norm_num in each stage or something like that?

[Mario Carneiro]

yes

[Kenny Lau]

it's definitely the wrong time for me to try, so maybe later

[Kenny Lau]

it's 2:28 AM

[Mario Carneiro]

I will come up with appropriate certificates in a bit, once I understand the situation right

[Floris van Doorn]

I'm currently adding a second argument and showing the monotonicity

[Floris van Doorn]

(updated the file on Github)

[Floris van Doorn]

are you working on the certificates @Mario Carneiro ? Otherwise I'll start on it now.

[Mario Carneiro]

I have this so far, for the lower bound:

{32 a5 >= 157/50, 2 - a5^2 >= a4, a4^2 >= 2 + a3, a3^2 >= 2 + a2,  a2^2 >= 2 + a1, a1^2 >= 2}

[Mario Carneiro]

Oh, I didn't expect this would work:

{a1 -> 2, a2 -> 2, a3 -> 2, a4 -> -(2687/1350), a5 -> 21/214}

[Mario Carneiro]

the last few (first few? my numbering is weird) are exactly 2, which works out nicely with the algebra. You probably want to special case that since it makes everything easy for going to higher values

[Floris van Doorn]

I'm still a bit confused what your variables a1-a5 correspond to.

[Mario Carneiro]

Here's the abstract argument:

pi >= 2^5 Sqrt[2 - f[0, 4]] >= 2^5 a5 >= 157/50
2 - a5^2 >= a4 >= f[a3, 1] >= f[a2, 2] >= f[a1, 3] >= f[0, 4]

[Mario Carneiro]

where f is your series function

[Mario Carneiro]

Note that there are no side conditions on this proof, because a4^2 >= 2 + a3 implies a4 >= Sqrt[2 + a3] regardless of whether a4 is positive

[Mario Carneiro]

for the upper bound I think there will be side conditions

[Floris van Doorn]

a4^2 >= 2 + a3 implies a4 >= Sqrt[2 + a3]

uhm... no: (-10) ^ 2 >= 2 + 0 but not -10 >= sqrt(2 + 0).

[Floris van Doorn]

For the upper bound there won't be side conditions, I indeed realized that a couple of minutes ago.

[Mario Carneiro]

Oh boy, I keep getting myself confused here. You are right, each of the a's must be positive, and I messed up a sign somewhere so my cert doesn't work

[Floris van Doorn]

The signs are pretty confusing indeed. But if you give me the certificate for the inequalities

f[a4, 0] >= f[a3, 1] >= f[a2, 2] >= f[a1, 3] >= f[0, 4]

then I'm ready to plug them in my Lean code now.

[Mario Carneiro]

Okay, here are the values for the lower bound:

{a1 -> 577/408, a2 -> 619/335, a3 -> 2144/1093, a4 -> 2687/1350, a5 -> 89/907}

[Floris van Doorn]

it works!

[Floris van Doorn]

lemma pi_gt_314 : pi > 3.14 :=
begin
  refine lt_of_le_of_lt _ (pi_gt_sqrt_two_add_series 4),
  rw [mul_comm],
  apply le_mul_of_div_le, norm_num, apply le_sqrt_of_sqr_le,
  rw [le_sub],
  refine le_trans (sqrt_two_add_series_step_up (577/408) (by norm_num) (by norm_num)) _,
  refine le_trans (sqrt_two_add_series_step_up (619/335) (by norm_num) (by norm_num)) _,
  refine le_trans (sqrt_two_add_series_step_up (2144/1093) (by norm_num) (by norm_num)) _,
  refine le_trans (sqrt_two_add_series_step_up (2687/1350) (by norm_num) (by norm_num)) _,
  rw [sqrt_two_add_series_zero],
  norm_num
end

[Floris van Doorn]

lemma pi_lt_315 : pi < 3.15 :=
begin
  refine lt_of_lt_of_le (pi_lt_sqrt_two_add_series 4) _,
  apply add_le_of_le_sub_right, rw [mul_comm], apply mul_le_of_le_div, apply pow_pos, norm_num,
  rw [sqrt_le_left, sub_le], swap, norm_num,
  refine le_trans _ (sqrt_two_add_series_step_down (41/29) (by norm_num)),
  refine le_trans _ (sqrt_two_add_series_step_down (194/105) (by norm_num)),
  refine le_trans _ (sqrt_two_add_series_step_down (13513/6889) (by norm_num)),
  refine le_trans _ (sqrt_two_add_series_step_down (2271/1141) (by norm_num)),
  norm_num
end

[Floris van Doorn]

The compiling proofs can now be found here: https://github.com/fpvandoorn/mathlib/blob/pi/src/data/real/pi.lean

They do each take ~5 seconds to compile. Maybe we can find smaller certificates (I found mine in a very naive way).

[Mario Carneiro]

I do want to figure out a general mechanism for optimizing all the certificates simultaneously; right now I'm optimizing them one by one and tweaking the results

[Johan Commelin]

Can't a tactic do that for you?

[Johan Commelin]

This must be a well-known idea in ITP, right?

[Mario Carneiro]

yes and no; this is stuff that should not be re-run every time

[Floris van Doorn]

Yes, I'm just using Mathematica's Rationalize myself, until Rationalize finds a rational number which is on the correct side of the number I'm approximating

[Johan Commelin]

The tactic could trace its result

[Mario Carneiro]

me too (you can ask for a number in a range with Rationalize[(a+b)/2,(a-b)/2])

[Floris van Doorn]

oh yes, good point

[Mario Carneiro]

there isn't much point doing this in lean, unless lean gets as good at numbers as a CAS

[Mario Carneiro]

this is really the sort of thing CASs are designed for: unverified computation / symbolic manipulation

[Mario Carneiro]

If I was to do this sort of thing a lot I would want to automate some parts of it

[Johan Commelin]

So we need a Sage bridge...

[Mario Carneiro]

sure, that would do it

[Mario Carneiro]

or code extraction

[Mario Carneiro]

in this case, it's just pure computation, although it's nice having mathematica with builtin functions like Rationalize so we don't have to write it on the spot

[Johan Commelin]

@Mario Carneiro Do you think meta code could at some point be just as fast a proper CAS? I don't see why we shouldn't strive for world domination there as well...

[Mario Carneiro]

sure, we just need more efficiency

[Johan Commelin]

Is there hope that :four_leaf_clover: will bring this efficiency?

[Mario Carneiro]

I've used lean plenty for this kind of thing, and it's a big bottleneck

[Mario Carneiro]

Yes, there is hope... I'm not putting all my eggs in that basket though

[Johan Commelin]

You haven't even put all your ITP eggs in Leans basket...

[Mario Carneiro]

I don't believe in world dominion

[Mario Carneiro]

it's not a healthy mindset. Better to support open communication and interaction between systems, because there will always be people that prefer another language

[Johan Commelin]

It also means you have to support n^2 communications...

[Mario Carneiro]

the common answer to that is "not if it's a star graph", i.e. world domination

[Mario Carneiro]

but other sparse graphs work too

[Mario Carneiro]

right now it's just an almost completely disconnected graph

[Kevin Buzzard]

I swear that was the happiest pi day of my life.

[Mario Carneiro]

this was a pretty good pi day activity

[Johan Commelin]

It would be cool to see some end result go into mathlib. Does that make sense?

[Mario Carneiro]

yes, I think it should certainly go in mathlib

[Scott Morrison]

Seems pretty reasonable. We can chuck out old estimates that have been superseded (both in terms of precision of the bounds and ideas), but the best available bounds sound like a useful thing.

[Mario Carneiro]

we might want to revisit the methods later, but something is better than nothing (and this method is pretty good already)

[Mario Carneiro]

I don't necessarily even think that we should only provide the best bounds. Tighter bounds are harder to calculate and harder to use

[Mario Carneiro]

In fact, as I argued in this thread, the best way to prove these kinds of inequalities is by controlled relaxation of the bounds to something simpler

[Kevin Buzzard]

https://xenaproject.wordpress.com/2019/03/15/happy-pi-day/

[Kevin Buzzard]

Comments welcome as ever.

[Johan Commelin]

@Kevin Buzzard Here is a link to the archive: https://leanprover-community.github.io/archive/113488general/89755happypiday.html
It doesn't require login

[Kenny Lau]

it's hard to generate good certs...

[Kenny Lau]

import math

# pi > 2 ^ (n+1) * sqrt (2 - sqrt_two_add_series 0 n)
# pi < 2 ^ (n+1) * sqrt (2 - sqrt_two_add_series 0 n) + 1 / 4 ^ n
# sqrt_two_add_series(x,n+1) = sqrt(2+sqrt_two_add_series(x,n))
def sqrt_two_add_series(x,n):
    while n: x=math.sqrt(2+x); n-=1
    return x

for n in range(10):
    print("n=%d:\t%.10f < pi < %.10f"%(n, 2**(n+1)*math.sqrt(2-sqrt_two_add_series(0,n)),
        2**(n+1)*math.sqrt(2-sqrt_two_add_series(0,n))+1/4**n))

"""OUTPUT:
n=0:    2.8284271247 < pi < 3.8284271247
n=1:    3.0614674589 < pi < 3.3114674589
n=2:    3.1214451523 < pi < 3.1839451523
n=3:    3.1365484905 < pi < 3.1521734905
n=4:    3.1403311570 < pi < 3.1442374070
n=5:    3.1412772509 < pi < 3.1422538134
n=6:    3.1415138011 < pi < 3.1417579418
n=7:    3.1415729404 < pi < 3.1416339755
n=8:    3.1415877253 < pi < 3.1416029841
n=9:    3.1415914215 < pi < 3.1415952362
"""

# to prove: 3.14 <=  2**5 * math.sqrt(2 - sqrt_two_add_series(0,4))
# to prove: 157/1600 <=  math.sqrt(2 - sqrt_two_add_series(0,4))
# to prove: 24649/2560000 <=  2 - sqrt_two_add_series(0,4)
# to prove: sqrt_two_add_series(0,4) <= 5095351/2560000
# rule: 2+x<=y^2 |- sqrt_two_add_series(x,n+1) <= sqrt_two_add_series(y,n)

print(sqrt_two_add_series(0,4))
print(sqrt_two_add_series(577/408,3))
print(sqrt_two_add_series(2537/1373,2))
print(sqrt_two_add_series(15163/7730,1))
print(sqrt_two_add_series(1447/727,0))
print(5095351/2560000)

"""OUTPUT:
1.9903694533443939
1.9903694901455364
1.9903707033912401
1.9903713890865138
1.9903713892709767
1.990371484375
"""

[Kevin Buzzard]

Use continued fractions

[Kevin Buzzard]

That is a cheap yet powerful method to approximate reals by rationals

[Kenny Lau]

that's what I used

[Kenny Lau]

here are some smaller certs:

print(sqrt_two_add_series(0,4))
print(sqrt_two_add_series(577/408,3))
print(sqrt_two_add_series(619/335,2))
print(sqrt_two_add_series(2297/1171,1))
print(sqrt_two_add_series(1447/727,0))
print(5095351/2560000)

"""OUTPUT:
1.9903694533443939
1.9903694901455364
1.9903695896706308
1.9903696406887783
1.9903713892709767
1.990371484375
"""

[Kenny Lau]

here are even smaller certs:

print(sqrt_two_add_series(0,4))
print(sqrt_two_add_series(577/408,3))
print(sqrt_two_add_series(619/335,2))
print(sqrt_two_add_series(1940/989,1))
print(sqrt_two_add_series(1447/727,0))
print(5095351/2560000)

"""OUTPUT:
1.9903694533443939
1.9903694901455364
1.9903695896706308
1.9903711590704518
1.9903713892709767
1.990371484375
"""

[Kenny Lau]

@Floris van Doorn I've made a PR to your file

[Kenny Lau]

even smaller:

print(sqrt_two_add_series(0,4))
print(sqrt_two_add_series(338/239,3))
print(sqrt_two_add_series(704/381,2))
print(sqrt_two_add_series(1940/989,1))
print(sqrt_two_add_series(1447/727,0))
print(2-(3.14/2**5)**2)

"""OUTPUT:
1.9903694533443939
1.990369667837741
1.9903700913568452
1.9903711590704518
1.9903713892709767
1.990371484375
"""

[Kenny Lau]

semi-brute-forced the (locally) smallest certs:

print(sqrt_two_add_series(0,4))
print(sqrt_two_add_series(99/70,3))
print(sqrt_two_add_series(874/473,2))
print(sqrt_two_add_series(1940/989,1))
print(sqrt_two_add_series(1447/727,0))
print(2-(3.14/2**5)**2)

"""OUTPUT:
1.9903694533443939
1.9903707035225453
1.9903708019895043
1.9903711590704518
1.9903713892709767
1.990371484375
"""

[Kenny Lau]

fitness = sum of all the numerators and denominators

[Kenny Lau]

def sqrt_two_add_series(x,n):
    while n: x=(2+x)**0.5; n-=1
    return x

def gcd(a,b):
    while a: a,b=b%a,a
    return b

def approx(lo,hi,N):
    res = []
    for denom in range(N):
        num = int(hi*denom)
        if gcd(num,denom) == 1 and lo <= num/denom:
            res.append((num,denom))
    return res

N=5000

res = [[L] for L in approx(sqrt_two_add_series(0,4), 2-(3.14/2**5)**2, N)]

for i in [3,2,1]:
    lo = sqrt_two_add_series(0,i)
    new = []
    for L in res:
        size = sum(num+denom for num,denom in L)
        hi = (L[0][0]/L[0][1])**2-2
        for num,denom in approx(lo,hi,N):
            if num+denom+size < 7000:
                new.append([(num,denom)] + L)
    res = new

min_size = 7000
for L in res:
    size = sum(num+denom for num,denom in L)
    if size < min_size:
        print(size, L)
        min_size = size

[Kenny Lau]

however the bottleneck is actually proving:

sqrt_two_add_series (1447 / 727) 0 ≤ 2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2

[Kenny Lau]

we didn't certify the right hand side

[Kevin Buzzard]

I think a Chinese mathematician 1000 years ago used 3•2^n-gons to estimate pi up to 7 digits?

@Kenny Lau do you think it's feasible to get up to 7 digits in Lean?

[Kenny Lau]

yes

[Kenny Lau]

@Mario Carneiro how would you prove ⊢ 1447 / 727 ≤ 2 - (157 / 50 / 2 ^ (4 + 1)) ^ 2?

[Kevin Buzzard]

Are these rationals or reals?

[Mario Carneiro]

There was a fifth cert in my version that floris omitted

[Kevin Buzzard]

I mean, I guess they're reals. Would it make life easier if they were rationals?

[Mario Carneiro]

you want to insert an intermediate: 1447 / 727 ≤ 2 - a^2 and a <= 157 / 50 / 2 ^ (4 + 1)

[Mario Carneiro]

good work on the cert optimization tho

[Kenny Lau]

I mean, I guess they're reals. Would it make life easier if they were rationals?

no, you almost never want to use the decidable instance

[Kenny Lau]

you want to insert an intermediate: 1447 / 727 ≤ 2 - a^2 and a <= 157 / 50 / 2 ^ (4 + 1)

ok lemme work on it

[Kevin Buzzard]

a=0 works great

[Mario Carneiro]

wait what? Something must be wrong then

[Mario Carneiro]

oh, a >= 157 / 50 / 2 ^ (4 + 1)

[Mario Carneiro]

like I mentioned before, sign headaches are real here

[Kevin Buzzard]

246/2507

[Kevin Buzzard]

@Kenny Lau

[Kenny Lau]

(157/50/2^5 = 157/1600 which works better, but actually this is not the intermediate we want)

[Kenny Lau]

lemma pi_gt_314 : pi > 3.14 :=
begin
  refine lt_of_le_of_lt _ (pi_gt_sqrt_two_add_series 4), rw [mul_comm],
  apply le_mul_of_div_le, norm_num, apply le_sqrt_of_sqr_le, rw [le_sub],
  refine le_trans (sqrt_two_add_series_step_up (99/70) (by norm_num) (by norm_num)) _,
  refine le_trans (sqrt_two_add_series_step_up (874/473) (by norm_num) (by norm_num)) _,
  refine le_trans (sqrt_two_add_series_step_up (1940/989) (by norm_num) (by norm_num)) _,
  refine le_trans (sqrt_two_add_series_step_up (1447/727) (by norm_num) (by norm_num)) _,
  refine le_trans (by norm_num : (1447 / 727 : ℝ) ≤ 2 - 141 / 14644) _,
  rw sub_le_sub_iff_left, norm_num
end

[Kenny Lau]

this is 2 seconds faster on my computer (6s vs 8s)

[Kenny Lau]

@Mario Carneiro do you see how to make this faster?

⊢ (157 / 50 / 2 ^ (4 + 1)) ^ 2 ≤ 141 / 14644

[Mario Carneiro]

same trick

[Mario Carneiro]

prove a^2 <= bla and bla <= a

[Mario Carneiro]

By the way, the actual optimization criterion is the number of theorem applications in the norm_num proofs, although this is crazy hard to predict

[Kenny Lau]

@Mario Carneiro but like bla would just be 141/14644

[Kenny Lau]

I already optimized this

[Mario Carneiro]

no, you are optimizing a

[Mario Carneiro]

157 / 50 / 2 ^ (4 + 1) <= magic and magic^2 <= 141/14644

[Kenny Lau]

the magic would just be 157/50/2^5 which is 157/1600

[Kenny Lau]

I already reached the limit I think

[Kenny Lau]

confirmed

[Kevin Buzzard]

(deleted)

[Mario Carneiro]

yeah okay then it's just a lot of working out

[Mario Carneiro]

If this is really expensive compared to the other steps, you could try making 141/14644 a bit larger and see if you can compensate with the others

[Kevin Buzzard]

Is clearing denominators expensive? There's a theorem which says you can do that

[Mario Carneiro]

Also I'm not sure but it might help to prepare the division cancellation steps since all of the magic numbers are fractions

[Mario Carneiro]

yeah that

[Kenny Lau]

what does that mean?

[Mario Carneiro]

meaning your intermediate theorems, the ones that you stick the magic numbers into, have side goals that talk only about multiplication and addition of natural numbers

[Kenny Lau]

what goal in particular?

[Mario Carneiro]

uh, I'm not sure what the exact form of the theorems being applied are

[Mario Carneiro]

but you can have the magic numbers be \u a / \u b instead of a real fraction

[Kenny Lau]

bottlenecks:

⊢ 2 + 0 ≤ (99 / 70) ^ 2
⊢ 2 + 99 / 70 ≤ (874 / 473) ^ 2
⊢ 2 + 874 / 473 ≤ (1940 / 989) ^ 2
⊢ 2 + 1940 / 989 ≤ (1447 / 727) ^ 2
⊢ 1447 / 727 ≤ 2 - 141 / 14644
⊢ (157 / 50 / 2 ^ (4 + 1)) ^ 2 ≤ 141 / 14644

[Mario Carneiro]

and then the side goals are things like a^2 * d <= (c + 2 d) * b^2 instead of (a/b)^2 <= c/d + 2

[Kenny Lau]

why would this be faster?

[Mario Carneiro]

Depending on how norm_num unfolds the things, it might work out larger intermediates than necessary

[Kenny Lau]

hmm

[Mario Carneiro]

In theory, norm_num should be doing the right thing, and so this won't actually gain much

[Mario Carneiro]

You might consider looking at the proof to profile it

[Mario Carneiro]

(split off the norm_num subproofs to lemmas)

[Kenny Lau]

the bottlenecks are essentially uncertifiable right?

[Mario Carneiro]

Well, at some point you actually have to prove a tight bound

[Kevin Buzzard]

Kenny how do you time how long a command takes?

[Mario Carneiro]

Generally speaking, the difficulty of the proof (size of the numbers) is inversely proportional to the size of the gap

[Mario Carneiro]

and we have placed 5 intermediate points in the original gap

[Mario Carneiro]

so if we shift the points around a bit we can equalize the load

[Kenny Lau]

for ⊢ 2 + 99 / 70 ≤ (874 / 473) ^ 2, the largest number norm_num produces in the proof is 361313348

[Mario Carneiro]

that sounds about right, for those numbers

[Kenny Lau]

(norm_num.subst_into_div (763876 / 473) 473 (874 * (874 / 473)) 473 (763876 / 223729)
                       (norm_num.div_eq_div_helper (763876 / 473) 473 763876 223729 361313348
                          (norm_num.subst_into_prod (763876 / 473) 223729 (763876 / 473) 223729 361313348
                             (norm_num.subst_into_div 763876 473 763876 473 (763876 / 473)
                                (norm_num.div_eq_div_helper 763876 473 763876 473 361313348
                                   (norm_num.subst_into_prod 763876 473 763876 473 361313348 (eq.refl 763876)
                                      (eq.refl 473)
                                      (norm_num.mul_bit1_helper 763876 236 180274736 361313348
                                         (norm_num.mul_bit0_helper 763876 118 90137368

[Mario Carneiro]

it should be 874^2 * 70

[Kenny Lau]

well it isn't

[Kevin Buzzard]

It's something times 473

[Kevin Buzzard]

it's 874^2 * 473

[Kenny Lau]

wat

[Kenny Lau]

how did that number get there

[Mario Carneiro]

I should note this is the one part of norm_num I didn't write

[Mario Carneiro]

I offload + - * / to leo's impl

[Mario Carneiro]

but it's not like there are that many ways to write a (correct) arithmetic prover

[Kenny Lau]

anyway I'm going now, this was fun

[Mario Carneiro]

good work

[Mario Carneiro]

aha:

variables {α : Type*} [field α]
#check @norm_num.div_eq_div_helper α _ (763876 / 473) 473 763876 223729 361313348
-- norm_num.div_eq_div_helper (763876 / 473) 473 763876 223729 361313348 :
--   763876 / 473 * 223729 = 361313348 →
--   763876 * 473 = 361313348 → 473 ≠ 0 → 223729 ≠ 0 → 763876 / 473 / 473 = 763876 / 223729

It is cross-multiplying to prove that (874 * (874 / 473)) / 473 = (763876 / 473) / 473 = 763876 / 223729

[Mario Carneiro]

that's... not a good idea for simplifying divisions of divisions

[Kevin Buzzard]

PR to core Kenny!

[Floris van Doorn]

At the bottom of the file https://github.com/fpvandoorn/mathlib/blob/pi/src/data/real/pi.lean I also computed pi up to 5 decimal digits. It does take 20-25 seconds to compile for me though.

[Floris van Doorn]

It's definitely possible to get better approximations, but the compilation time will just be bigger for that. So I don't think we want the tightest bounds in mathlib. But soon I'll PR the first 3 digits of pi to mathlib.

[Mario Carneiro]

what is the profiler result? Is it norm_num or the kernel taking up the time?

[Floris van Doorn]

Not sure how to read the output, here it is for pi_gt_31415:

elaboration of pi_gt_31415 took 11.9s

elaboration: tactic execution took 8.56s
num. allocated objects:  8023
num. allocated closures: 9813
num. allocated big nums: 96
 8560ms   100.0%   tactic.step
 8560ms   100.0%   scope_trace
 8560ms   100.0%   tactic.istep
 8560ms   100.0%   _interaction._lambda_2
 8560ms   100.0%   tactic.istep._lambda_1
 5972ms    69.8%   tactic.refine
 5945ms    69.5%   tactic.to_expr
 2583ms    30.2%   interaction_monad.orelse'
 2580ms    30.1%   tactic.seq
 1641ms    19.2%   interaction_monad_orelse
 1621ms    18.9%   tactic.interactive.norm_num1
 1620ms    18.9%   tactic.replace_at._lambda_4
 1620ms    18.9%   tactic.alternative._lambda_3
 1620ms    18.9%   tactic.replace_at
 1615ms    18.9%   interaction_monad.monad._lambda_9
 1614ms    18.9%   tactic.ext_simplify_core
 1614ms    18.9%   norm_num.derive
 1607ms    18.8%   tactic.repeat
 1607ms    18.8%   _private.3096317357.repeat_aux._main._lambda_1
 1607ms    18.8%   tactic.all_goals
 1607ms    18.8%   _private.3321893233.all_goals_core._main._lambda_2
 1607ms    18.8%   tactic.try_core
 1607ms    18.8%   all_goals_core
 1607ms    18.8%   repeat_aux
 1548ms    18.1%   norm_num.derive._main._lambda_3
 1364ms    15.9%   tactic.norm_num
  959ms    11.2%   tactic.interactive.simp_core
  937ms    10.9%   tactic.join_user_simp_lemmas
  937ms    10.9%   tactic.mk_simp_set_core
  937ms    10.9%   simp_lemmas.mk_default
  937ms    10.9%   tactic.mk_simp_set
[...]

[Floris van Doorn]

So most of the tactic time is in refine, so building the proof term(?). But there is also 3-4 seconds of time which is not used to execute tactics. Is that just to type-check the final proof?

[Mario Carneiro]

I think you missed out on some part... there should be a line for kernel typechecking

[Mario Carneiro]

it's next to the elaboration of foo took n s line

[Floris van Doorn]

Yeah, I thought so too. I cannot find it anywhere in my Lean messages

[Mario Carneiro]

maybe that just means it was below the threshold

[Mario Carneiro]

but damn that's a lot of to_expr

[Mario Carneiro]

and istep too... the tactic overhead is significant

[Floris van Doorn]

Oh no, I know what is going on: all the tactic.refine are the nested norm_nums

[Mario Carneiro]

The way you read this is to look for nodes that come just before a percentage drop

[Floris van Doorn]

When we execute

refine sqrt_two_add_series_step_up (11482/8119)  (by norm_num) (by norm_num) _,

the inner norm_num reports that the inner norm_num takes up almost a second:

713ms    98.5%   norm_num.derive

but for the outer tactic block, this is all part of the refine tactic.

[Mario Carneiro]

refine is quick, it is just calling another tactic that is slow

[Floris van Doorn]

ok, sure. But the to_expr just calls the nested norm_num tactics, right?

[Mario Carneiro]

yes, but that should show up in the total profile

[Mario Carneiro]

1614ms    18.9%   norm_num.derive
1607ms    18.8%   tactic.repeat
...

[Mario Carneiro]

compared to 8560ms for all the tactics, that's not as large a fraction as I would have hoped

[Floris van Doorn]

Mario Carneiro: yes, but that should show up in the total profile

I don't think it does. The sum of the times spend in norm_num.derive in the nested by norm_num tactics is much larger than the reported time on norm_num.derive when I hover over the end of the tactic block

[Floris van Doorn]

If I add everything up, about 5.3s is spend in norm_num.derive

[Floris van Doorn]

A bit more actually, but between 5s and 6s.

[Mario Carneiro]

It's surprising how inconclusive this test is, but it looks like the outer block does in fact reflect slow downs caused by inner blocks

set_option profiler true

example : true :=
begin
  tactic.sleep 10000,
  refine begin
    tactic.sleep 10000,
    trivial
  end
end

[Mario Carneiro]

either the profiler is way off, or sleep is not very accurate

[Mario Carneiro]

elaboration: tactic compilation took 17.8ms
test.lean:10:0: information tactic profile data
elaboration: tactic execution took 14.4s
num. allocated objects:  15
num. allocated closures: 46
14404ms   100.0%   tactic.istep
14404ms   100.0%   _interaction._lambda_2
14404ms   100.0%   scope_trace
14404ms   100.0%   tactic.istep._lambda_1
14404ms   100.0%   tactic.step
 7470ms    51.9%   tactic.to_expr
 7470ms    51.9%   tactic.refine
 6934ms    48.1%   tactic.sleep

[Floris van Doorn]

That seems pretty conclusive to me... The bottom end reports 1.7s spend at sleep and 5.2s at refine, the top end reports 5.2s spend at sleep. But yeah, those numbers are not very close to 2000 and 4000...

[Floris van Doorn]

(for future reference: the sleep numbers in Mario's post were originally 2000 and 4000)

[Mario Carneiro]

even when I crank the numbers up they are still not so close

[Mario Carneiro]

I think the units must be off

[Floris van Doorn]

yeah, there must be something off with sleep. But it's clear that the inner sleep time only counts towards refine, not towards sleep.

[Mario Carneiro]

yeah. To make it even more obvious:

example :=
begin
  have := begin
    tactic.sleep 2000,
    trivial
  end,
  refine begin
    tactic.sleep 2000,
    trivial
  end
end

elaboration: tactic execution took 2.8s
num. allocated objects:  21
num. allocated closures: 52
 2798ms   100.0%   tactic.to_expr
 2798ms   100.0%   tactic.istep
 2798ms   100.0%   _interaction._lambda_2
 2798ms   100.0%   scope_trace
 2798ms   100.0%   tactic.istep._lambda_1
 2798ms   100.0%   tactic.step
 1429ms    51.1%   tactic.interactive.have._lambda_1
 1369ms    48.9%   tactic.refine

[Mario Carneiro]

There is an easy solution to this for your situation though: don't use nested tactics, use a refine with more underscores and kill the subgoals with norm_num in the main tactic block

[Floris van Doorn]

If I do that, we see that indeed most of the time is spend at norm_num:

elaboration: tactic execution took 7.48s
num. allocated objects:  29602
num. allocated closures: 30223
num. allocated big nums: 322
 7482ms   100.0%   tactic.istep._lambda_1
 7482ms   100.0%   _interaction._lambda_2
 7482ms   100.0%   tactic.istep
 7482ms   100.0%   tactic.step
 7482ms   100.0%   scope_trace
 7442ms    99.5%   interaction_monad.orelse'
 7439ms    99.4%   tactic.seq
 7420ms    99.2%   interaction_monad.monad._lambda_9
 7388ms    98.7%   tactic.all_goals
 7388ms    98.7%   _private.3321893233.all_goals_core._main._lambda_2
 7388ms    98.7%   all_goals_core
 6640ms    88.7%   interaction_monad_orelse
 6619ms    88.5%   tactic.replace_at._lambda_4
 6619ms    88.5%   tactic.replace_at
 6619ms    88.5%   tactic.alternative._lambda_3
 6619ms    88.5%   tactic.interactive.norm_num1
 6563ms    87.7%   norm_num.derive
 6563ms    87.7%   tactic.ext_simplify_core
 6499ms    86.9%   norm_num.derive._main._lambda_3
 5696ms    76.1%   tactic.norm_num
 2613ms    34.9%   norm_num.eval_ineq
 2608ms    34.9%   norm_num.prove_lt
 1357ms    18.1%   _private.3096317357.repeat_aux._main._lambda_1
[...]

[Floris van Doorn]

Interestingly, there is a drop after all_goals_core (I'm now accumulating all the inequalities as goals, and killing them all with all_goals {norm_num} at the end)

[Floris van Doorn]

oh, in the other proof (pi_lt_31416) the drop-offs from 100% to 90% are at very different tactics. So maybe that is just noise/semi-randomness.

[Floris van Doorn]

Since people requested it, I also computed the first 7 digits of pi at the bottom of https://github.com/fpvandoorn/mathlib/blob/pi/src/data/real/pi.lean