## Stream: general

### Topic: happy pi day!

happy pi day!

#### Kevin Buzzard (Mar 14 2019 at 12:41):

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

#### Kevin Buzzard (Mar 14 2019 at 12:54):

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 (Mar 14 2019 at 13:27):

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 (Mar 14 2019 at 13:30):

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 (Mar 14 2019 at 13:32):

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 (Mar 14 2019 at 13:38):

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 (Mar 14 2019 at 13:40):

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

#### Kenny Lau (Mar 14 2019 at 13:40):

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

#### Chris Hughes (Mar 14 2019 at 14:04):

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 (Mar 14 2019 at 14:23):

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 (Mar 14 2019 at 14:26):

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 (Mar 14 2019 at 14:43):

that power series is notoriously slow to converge

#### Kevin Buzzard (Mar 14 2019 at 14:44):

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 (Mar 14 2019 at 14:45):

Which is more than we know currently in mathlib.

#### Kevin Buzzard (Mar 14 2019 at 14:47):

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 (Mar 14 2019 at 14:47):

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

#### Kevin Buzzard (Mar 14 2019 at 14:48):

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

#### Kevin Buzzard (Mar 14 2019 at 14:48):

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 (Mar 14 2019 at 14:48):

You'll need pi as an integral, correct?

#### Kevin Buzzard (Mar 14 2019 at 14:49):

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 (Mar 14 2019 at 14:49):

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

#### Abhimanyu Pallavi Sudhir (Mar 14 2019 at 14:50):

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

ha ha

#### Kevin Buzzard (Mar 14 2019 at 14:50):

there's a problem with the polygon approach

#### Kevin Buzzard (Mar 14 2019 at 14:51):

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

#### Kevin Buzzard (Mar 14 2019 at 14:52):

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

#### Kevin Buzzard (Mar 14 2019 at 14:52):

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

#### Kevin Buzzard (Mar 14 2019 at 14:53):

and that the derivative of sin was cos

#### Kevin Buzzard (Mar 14 2019 at 14:53):

then setting n=6 you get pi>3

#### Reid Barton (Mar 14 2019 at 14:53):

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

#### Kevin Buzzard (Mar 14 2019 at 14:54):

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

#### Kevin Buzzard (Mar 14 2019 at 14:54):

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 (Mar 14 2019 at 14:54):

that crazy quote bug again

#### Kevin Buzzard (Mar 14 2019 at 14:55):

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

#### Kevin Buzzard (Mar 14 2019 at 14:55):

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 (Mar 14 2019 at 14:56):

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

#### Kevin Buzzard (Mar 14 2019 at 14:56):

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

#### Kevin Buzzard (Mar 14 2019 at 14:56):

not just some blather about cloud computing

#### Kevin Buzzard (Mar 14 2019 at 14:57):

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

#### Kenny Lau (Mar 14 2019 at 15:01):

firstly how do you define pi

it's in mathlib

#### Kevin Buzzard (Mar 14 2019 at 15:02):

so it is defined by official decree already

#### Kevin Buzzard (Mar 14 2019 at 15:02):

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

it's in mathlib?

#### Kevin Buzzard (Mar 14 2019 at 15:03):

no wait, that was something else.

#### Kevin Buzzard (Mar 14 2019 at 15:03):

Isn't it in mathlib??

yes it's there

#### Kevin Buzzard (Mar 14 2019 at 15:04):

It's in analysis.exponential

#### Kevin Buzzard (Mar 14 2019 at 15:04):

local notation π := 22 / 7

no wait

#### Kenny Lau (Mar 14 2019 at 15:05):

oh right they proved sin(pi) = 0 right

#### Kevin Buzzard (Mar 14 2019 at 15:05):

local notation π := real.pi

oh

#### Kenny Lau (Mar 14 2019 at 15:05):

what's the definition?

#### Abhimanyu Pallavi Sudhir (Mar 14 2019 at 15:05):

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 (Mar 14 2019 at 15:06):

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

#### Kevin Buzzard (Mar 14 2019 at 15:07):

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 (Mar 14 2019 at 15:07):

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

#### Kevin Buzzard (Mar 14 2019 at 15:08):

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

#### Kevin Buzzard (Mar 14 2019 at 15:08):

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

#### Kevin Buzzard (Mar 14 2019 at 15:08):

Maybe it's not less than 4 after all!

#### Kevin Buzzard (Mar 14 2019 at 15:08):

pesky axiom of choice

#### Kevin Buzzard (Mar 14 2019 at 15:10):

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


#### Kevin Buzzard (Mar 14 2019 at 15:10):

Proof uses cos_one_pos and cos_two_neg and IVT

#### Kevin Buzzard (Mar 14 2019 at 15:10):

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

#### Kevin Buzzard (Mar 14 2019 at 15:11):

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 (Mar 14 2019 at 15:12):

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 (Mar 14 2019 at 15:12):

Do we know the derivative of exp is exp yet?

#### Abhimanyu Pallavi Sudhir (Mar 14 2019 at 15:14):

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

#### Reid Barton (Mar 14 2019 at 15:14):

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 (Mar 14 2019 at 15:15):

cos_lt_cos_of_nonneg_of_le_pi_div_two

#### Reid Barton (Mar 14 2019 at 15:15):

or actually cos_lt_cos_of_nonneg_of_le_pi

#### Reid Barton (Mar 14 2019 at 15:16):

It's proved using angle addition formulas

#### Floris van Doorn (Mar 14 2019 at 15:49):

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.

nice!

#### Chris Hughes (Mar 14 2019 at 15:57):

I can imagine this being a long annoying proof.

#### Reid Barton (Mar 14 2019 at 15:57):

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

#### Floris van Doorn (Mar 14 2019 at 15:58):

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

Happy pi day :-)

#### Kevin Buzzard (Mar 14 2019 at 16:03):

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 (Mar 14 2019 at 16:03):

Or is there something else you're worried about?

#### Floris van Doorn (Mar 14 2019 at 16:03):

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

#### Kevin Buzzard (Mar 14 2019 at 16:03):

These things should go in pi.lean

#### Kevin Buzzard (Mar 14 2019 at 16:03):

norm_num proves 2 ne 0

#### Kevin Buzzard (Mar 14 2019 at 16:03):

or if it doesn't then take real parts

#### Kevin Buzzard (Mar 14 2019 at 16:03):

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

#### Chris Hughes (Mar 14 2019 at 16:03):

two_ne_zero'

rofl

#### Chris Hughes (Mar 14 2019 at 16:04):

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

#### Floris van Doorn (Mar 14 2019 at 16:04):

yes, I found the other one. Thanks

#### Kevin Buzzard (Mar 14 2019 at 16:04):

linear ordered -> char 0

#### Chris Hughes (Mar 14 2019 at 16:05):

The linear order version is in core I think.

#### Kevin Buzzard (Mar 14 2019 at 16:07):

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

#### Kevin Buzzard (Mar 14 2019 at 16:07):

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

#### Floris van Doorn (Mar 14 2019 at 19:06):

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 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,
norm_num, norm_num, norm_num,
rw [sub_nonneg, sqrt_le_left], norm_num, norm_num
end

end real


#### Kevin Buzzard (Mar 14 2019 at 19:08):

This doesn't compile for me

#### Kevin Buzzard (Mar 14 2019 at 19:08):

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


We need pi <= 8

#### Kevin Buzzard (Mar 14 2019 at 19:09):

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,
norm_num, norm_num, norm_num,
rw [sub_nonneg, sqrt_le_left], norm_num, norm_num
end


#### Floris van Doorn (Mar 14 2019 at 19:09):

oops, I missed that. One minute

#### Floris van Doorn (Mar 14 2019 at 19:11):

I updated the previous post, it should compile now.

#### Reid Barton (Mar 14 2019 at 19:11):

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 (Mar 14 2019 at 19:14):

@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 (Mar 14 2019 at 19:15):

I updated the code once more.

#### Kevin Buzzard (Mar 14 2019 at 19:29):

This should all go in pi.lean, right?

#### Kevin Buzzard (Mar 14 2019 at 19:30):

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

#### Kevin Buzzard (Mar 14 2019 at 19:30):

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

#### Kevin Buzzard (Mar 14 2019 at 19:31):

ie floor(pi*10)%10=1

#### Johan Commelin (Mar 14 2019 at 19:31):

This should all go in pi.lean, right?

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

#### Floris van Doorn (Mar 14 2019 at 19:37):

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 (Mar 14 2019 at 19:45):

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

probably not

#### Floris van Doorn (Mar 14 2019 at 23:35):

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 (Mar 14 2019 at 23:45):

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 (Mar 14 2019 at 23:57):

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

#### Kevin Buzzard (Mar 14 2019 at 23:57):

It's not even close

#### Kevin Buzzard (Mar 14 2019 at 23:58):

but I'm struggling with the last bit.

#### Kevin Buzzard (Mar 14 2019 at 23:58):

oh -- these are reals

#### Kevin Buzzard (Mar 14 2019 at 23:59):

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 (Mar 15 2019 at 00:05):

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 (Mar 15 2019 at 00:08):

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 (Mar 15 2019 at 00:08):

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

#### Kevin Buzzard (Mar 15 2019 at 00:20):

Only 31999999999998 digits to go and we have the record!

oh boy

#### Floris van Doorn (Mar 15 2019 at 01:49):

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 (Mar 15 2019 at 01:54):

I've verified all of your inequalities on WolframAlpha

#### Mario Carneiro (Mar 15 2019 at 01:56):

Oh boy I've missed a lot

#### Mario Carneiro (Mar 15 2019 at 01:57):

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

#### Mario Carneiro (Mar 15 2019 at 01:58):

Are all the inequalities implied by one of them?

#### Kenny Lau (Mar 15 2019 at 01:58):

the first one is exceptionally tight

#### Kenny Lau (Mar 15 2019 at 01:58):

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

#### Kenny Lau (Mar 15 2019 at 01:58):

wait are those approximations to sqrt(2)?

#### Kenny Lau (Mar 15 2019 at 01:59):

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

#### Kenny Lau (Mar 15 2019 at 01:59):

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

#### Mario Carneiro (Mar 15 2019 at 01:59):

I think they are approximations to 2

#### Mario Carneiro (Mar 15 2019 at 01:59):

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

interesting...?

#### Mario Carneiro (Mar 15 2019 at 02:01):

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

   1128227574827648531222203602074520069222725565123913061124308255143227446400872965401873904861225601 /
3402823669209384634633746074317682114560000000000000000000000000000000000000000000000000000000000000000
>=0


#### Mario Carneiro (Mar 15 2019 at 02:01):

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

#### Mario Carneiro (Mar 15 2019 at 02:02):

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

aha

#### Mario Carneiro (Mar 15 2019 at 02:05):

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

Where did this come from?

#### Floris van Doorn (Mar 15 2019 at 02:06):

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

#### Kenny Lau (Mar 15 2019 at 02:08):

/- 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

lemma pi_lt_315 : pi < 3.15 :=
begin


extract ^

#### Kenny Lau (Mar 15 2019 at 02:09):

anyway @Floris van Doorn remarkable work!

#### Mario Carneiro (Mar 15 2019 at 02:11):

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

#### Floris van Doorn (Mar 15 2019 at 02:12):

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 (Mar 15 2019 at 02:12):

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 (Mar 15 2019 at 02:14):

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

#### Floris van Doorn (Mar 15 2019 at 02:14):

The defining equation of sqrt_two_add_series is

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


#### Mario Carneiro (Mar 15 2019 at 02:14):

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

#### Floris van Doorn (Mar 15 2019 at 02:15):

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

#### Kenny Lau (Mar 15 2019 at 02:16):

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

#### Mario Carneiro (Mar 15 2019 at 02:17):

blasty lean users aren't going to like the answer

#### Mario Carneiro (Mar 15 2019 at 02:17):

you give intermediate small rational numbers as certificates

#### Floris van Doorn (Mar 15 2019 at 02:18):

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 (Mar 15 2019 at 02:19):

fortunately I'm not a blasty lean user

#### Mario Carneiro (Mar 15 2019 at 02:19):

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 (Mar 15 2019 at 02:19):

sorry playing catchup

#### Floris van Doorn (Mar 15 2019 at 02:21):

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

#### Mario Carneiro (Mar 15 2019 at 02:21):

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 (Mar 15 2019 at 02:22):

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 (Mar 15 2019 at 02:23):

and those numbers are the certificate

that sounds good

#### Floris van Doorn (Mar 15 2019 at 02:23):

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

#### Mario Carneiro (Mar 15 2019 at 02:24):

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

#### Mario Carneiro (Mar 15 2019 at 02:24):

and pick n so it works

#### Mario Carneiro (Mar 15 2019 at 02:25):

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

#### Mario Carneiro (Mar 15 2019 at 02:26):

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

#### Kenny Lau (Mar 15 2019 at 02:26):

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

#### Mario Carneiro (Mar 15 2019 at 02:26):

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

#### Kenny Lau (Mar 15 2019 at 02:26):

would that be useful?

#### Mario Carneiro (Mar 15 2019 at 02:27):

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

#### Kenny Lau (Mar 15 2019 at 02:28):

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

yes

#### Kenny Lau (Mar 15 2019 at 02:28):

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

it's 2:28 AM

#### Mario Carneiro (Mar 15 2019 at 02:29):

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

#### Floris van Doorn (Mar 15 2019 at 02:40):

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

#### Floris van Doorn (Mar 15 2019 at 02:51):

(updated the file on Github)

#### Floris van Doorn (Mar 15 2019 at 03:17):

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

#### Mario Carneiro (Mar 15 2019 at 03:18):

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 (Mar 15 2019 at 03:26):

Oh, I didn't expect this would work:

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


#### Mario Carneiro (Mar 15 2019 at 03:27):

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 (Mar 15 2019 at 03:27):

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

#### Mario Carneiro (Mar 15 2019 at 03:28):

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 (Mar 15 2019 at 03:28):

where f is your series function

#### Mario Carneiro (Mar 15 2019 at 03:31):

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 (Mar 15 2019 at 03:31):

for the upper bound I think there will be side conditions

#### Floris van Doorn (Mar 15 2019 at 03:37):

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

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

#### Floris van Doorn (Mar 15 2019 at 03:37):

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

#### Mario Carneiro (Mar 15 2019 at 03:45):

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 (Mar 15 2019 at 03:46):

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 (Mar 15 2019 at 04:02):

Okay, here are the values for the lower bound:

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


it works!

#### Floris van Doorn (Mar 15 2019 at 04:07):

lemma pi_gt_314 : pi > 3.14 :=
begin
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)) _,
norm_num
end


#### Floris van Doorn (Mar 15 2019 at 04:19):

lemma pi_lt_315 : pi < 3.15 :=
begin
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 (Mar 15 2019 at 04:22):

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 (Mar 15 2019 at 04:23):

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 (Mar 15 2019 at 04:23):

Can't a tactic do that for you?

#### Johan Commelin (Mar 15 2019 at 04:23):

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

#### Mario Carneiro (Mar 15 2019 at 04:23):

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

#### Floris van Doorn (Mar 15 2019 at 04:23):

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 (Mar 15 2019 at 04:24):

The tactic could trace its result

#### Mario Carneiro (Mar 15 2019 at 04:24):

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

#### Floris van Doorn (Mar 15 2019 at 04:25):

oh yes, good point

#### Mario Carneiro (Mar 15 2019 at 04:26):

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

#### Mario Carneiro (Mar 15 2019 at 04:27):

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

#### Mario Carneiro (Mar 15 2019 at 04:27):

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

#### Johan Commelin (Mar 15 2019 at 04:34):

So we need a Sage bridge...

#### Mario Carneiro (Mar 15 2019 at 04:38):

sure, that would do it

#### Mario Carneiro (Mar 15 2019 at 04:39):

or code extraction

#### Mario Carneiro (Mar 15 2019 at 04:39):

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 (Mar 15 2019 at 04:52):

@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 (Mar 15 2019 at 04:52):

sure, we just need more efficiency

#### Johan Commelin (Mar 15 2019 at 04:52):

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

#### Mario Carneiro (Mar 15 2019 at 04:52):

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

#### Mario Carneiro (Mar 15 2019 at 04:53):

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

#### Mario Carneiro (Mar 15 2019 at 04:54):

I don't believe in world dominion

#### Mario Carneiro (Mar 15 2019 at 04:55):

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 (Mar 15 2019 at 04:56):

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

#### Mario Carneiro (Mar 15 2019 at 04:57):

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

#### Mario Carneiro (Mar 15 2019 at 04:57):

but other sparse graphs work too

#### Mario Carneiro (Mar 15 2019 at 04:57):

right now it's just an almost completely disconnected graph

#### Kevin Buzzard (Mar 15 2019 at 08:39):

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

#### Mario Carneiro (Mar 15 2019 at 08:42):

this was a pretty good pi day activity

#### Johan Commelin (Mar 15 2019 at 08:57):

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

#### Mario Carneiro (Mar 15 2019 at 09:02):

yes, I think it should certainly go in mathlib

#### Scott Morrison (Mar 15 2019 at 09:02):

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 (Mar 15 2019 at 09:03):

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

#### Mario Carneiro (Mar 15 2019 at 09:03):

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 (Mar 15 2019 at 09:04):

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 (Mar 15 2019 at 09:31):

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

#### Johan Commelin (Mar 15 2019 at 09:51):

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

#### Kenny Lau (Mar 15 2019 at 09:53):

it's hard to generate good certs...

#### Kenny Lau (Mar 15 2019 at 09:53):

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
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)),

"""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

print(5095351/2560000)

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


#### Kevin Buzzard (Mar 15 2019 at 09:58):

Use continued fractions

#### Kevin Buzzard (Mar 15 2019 at 09:58):

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

#### Kenny Lau (Mar 15 2019 at 10:01):

that's what I used

#### Kenny Lau (Mar 15 2019 at 10:03):

here are some smaller certs:

print(sqrt_two_add_series(0,4))
print(5095351/2560000)

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


#### Kenny Lau (Mar 15 2019 at 10:08):

here are even smaller certs:

print(sqrt_two_add_series(0,4))
print(5095351/2560000)

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


#### Kenny Lau (Mar 15 2019 at 10:48):

even smaller:

print(sqrt_two_add_series(0,4))
print(2-(3.14/2**5)**2)

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


#### Kenny Lau (Mar 15 2019 at 11:05):

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

print(sqrt_two_add_series(0,4))
print(2-(3.14/2**5)**2)

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


#### Kenny Lau (Mar 15 2019 at 11:05):

fitness = sum of all the numerators and denominators

#### Kenny Lau (Mar 15 2019 at 11:09):

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]:
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 (Mar 15 2019 at 11:12):

however the bottleneck is actually proving:

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


#### Kenny Lau (Mar 15 2019 at 11:12):

we didn't certify the right hand side

#### Kevin Buzzard (Mar 15 2019 at 11:13):

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?

yes

#### Kenny Lau (Mar 15 2019 at 11:22):

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

#### Kevin Buzzard (Mar 15 2019 at 11:24):

Are these rationals or reals?

#### Mario Carneiro (Mar 15 2019 at 11:24):

There was a fifth cert in my version that floris omitted

#### Kevin Buzzard (Mar 15 2019 at 11:24):

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

#### Mario Carneiro (Mar 15 2019 at 11:25):

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

#### Mario Carneiro (Mar 15 2019 at 11:26):

good work on the cert optimization tho

#### Kenny Lau (Mar 15 2019 at 11:26):

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 (Mar 15 2019 at 11:26):

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

ok lemme work on it

a=0 works great

#### Mario Carneiro (Mar 15 2019 at 11:29):

wait what? Something must be wrong then

#### Mario Carneiro (Mar 15 2019 at 11:30):

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

#### Mario Carneiro (Mar 15 2019 at 11:30):

like I mentioned before, sign headaches are real here

#### Kevin Buzzard (Mar 15 2019 at 11:31):

246/2507

@Kenny Lau

#### Kenny Lau (Mar 15 2019 at 11:33):

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

#### Kenny Lau (Mar 15 2019 at 11:33):

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 (Mar 15 2019 at 11:34):

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

#### Kenny Lau (Mar 15 2019 at 11:35):

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

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


same trick

#### Mario Carneiro (Mar 15 2019 at 11:35):

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

#### Mario Carneiro (Mar 15 2019 at 11:37):

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 (Mar 15 2019 at 11:39):

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

#### Mario Carneiro (Mar 15 2019 at 11:39):

no, you are optimizing a

#### Mario Carneiro (Mar 15 2019 at 11:40):

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

#### Kenny Lau (Mar 15 2019 at 11:40):

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

#### Kenny Lau (Mar 15 2019 at 11:40):

I already reached the limit I think

confirmed

(deleted)

#### Mario Carneiro (Mar 15 2019 at 11:41):

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

#### Mario Carneiro (Mar 15 2019 at 11:42):

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 (Mar 15 2019 at 11:44):

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

#### Mario Carneiro (Mar 15 2019 at 11:44):

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

yeah that

#### Kenny Lau (Mar 15 2019 at 11:44):

what does that mean?

#### Mario Carneiro (Mar 15 2019 at 11:45):

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 (Mar 15 2019 at 11:46):

what goal in particular?

#### Mario Carneiro (Mar 15 2019 at 11:47):

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

#### Mario Carneiro (Mar 15 2019 at 11:47):

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

#### Kenny Lau (Mar 15 2019 at 11:49):

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 (Mar 15 2019 at 11:49):

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 (Mar 15 2019 at 11:49):

why would this be faster?

#### Mario Carneiro (Mar 15 2019 at 11:50):

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

hmm

#### Mario Carneiro (Mar 15 2019 at 11:51):

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

#### Mario Carneiro (Mar 15 2019 at 11:51):

You might consider looking at the proof to profile it

#### Mario Carneiro (Mar 15 2019 at 11:52):

(split off the norm_num subproofs to lemmas)

#### Kenny Lau (Mar 15 2019 at 11:54):

the bottlenecks are essentially uncertifiable right?

#### Mario Carneiro (Mar 15 2019 at 11:54):

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

#### Kevin Buzzard (Mar 15 2019 at 11:54):

Kenny how do you time how long a command takes?

#### Mario Carneiro (Mar 15 2019 at 11:55):

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

#### Mario Carneiro (Mar 15 2019 at 11:55):

and we have placed 5 intermediate points in the original gap

#### Mario Carneiro (Mar 15 2019 at 11:55):

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

#### Kenny Lau (Mar 15 2019 at 11:55):

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

#### Mario Carneiro (Mar 15 2019 at 11:56):

that sounds about right, for those numbers

#### Kenny Lau (Mar 15 2019 at 11:56):

(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 (Mar 15 2019 at 11:56):

it should be 874^2 * 70

well it isn't

#### Kevin Buzzard (Mar 15 2019 at 11:56):

It's something times 473

it's 874^2 * 473

wat

#### Kenny Lau (Mar 15 2019 at 11:57):

how did that number get there

#### Mario Carneiro (Mar 15 2019 at 11:58):

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

#### Mario Carneiro (Mar 15 2019 at 11:58):

I offload + - * / to leo's impl

#### Mario Carneiro (Mar 15 2019 at 11:59):

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

#### Kenny Lau (Mar 15 2019 at 12:00):

anyway I'm going now, this was fun

good work

#### Mario Carneiro (Mar 15 2019 at 12:04):

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 (Mar 15 2019 at 12:04):

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

#### Kevin Buzzard (Mar 15 2019 at 12:06):

PR to core Kenny!

#### Floris van Doorn (Mar 15 2019 at 13:38):

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 (Mar 15 2019 at 13:39):

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 (Mar 15 2019 at 14:16):

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

#### Floris van Doorn (Mar 15 2019 at 14:22):

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
2580ms    30.1%   tactic.seq
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
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 (Mar 15 2019 at 14:24):

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 (Mar 15 2019 at 14:30):

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

#### Mario Carneiro (Mar 15 2019 at 14:30):

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

#### Floris van Doorn (Mar 15 2019 at 14:30):

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

#### Mario Carneiro (Mar 15 2019 at 14:31):

maybe that just means it was below the threshold

#### Mario Carneiro (Mar 15 2019 at 14:31):

but damn that's a lot of to_expr

#### Mario Carneiro (Mar 15 2019 at 14:32):

and istep too... the tactic overhead is significant

#### Floris van Doorn (Mar 15 2019 at 14:32):

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

#### Mario Carneiro (Mar 15 2019 at 14:33):

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

#### Floris van Doorn (Mar 15 2019 at 14:33):

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 (Mar 15 2019 at 14:34):

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

#### Floris van Doorn (Mar 15 2019 at 14:34):

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

#### Mario Carneiro (Mar 15 2019 at 14:35):

yes, but that should show up in the total profile

#### Mario Carneiro (Mar 15 2019 at 14:35):

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


#### Mario Carneiro (Mar 15 2019 at 14:36):

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

#### Floris van Doorn (Mar 15 2019 at 14:37):

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 (Mar 15 2019 at 14:38):

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

#### Floris van Doorn (Mar 15 2019 at 14:39):

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

#### Mario Carneiro (Mar 15 2019 at 14:44):

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 (Mar 15 2019 at 14:44):

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

#### Mario Carneiro (Mar 15 2019 at 14:46):

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 (Mar 15 2019 at 14:46):

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 (Mar 15 2019 at 14:47):

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

#### Mario Carneiro (Mar 15 2019 at 14:47):

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

#### Mario Carneiro (Mar 15 2019 at 14:47):

I think the units must be off

#### Floris van Doorn (Mar 15 2019 at 14:48):

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 (Mar 15 2019 at 14:49):

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 (Mar 15 2019 at 14:52):

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 (Mar 15 2019 at 15:00):

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
7439ms    99.4%   tactic.seq
7388ms    98.7%   tactic.all_goals
7388ms    98.7%   _private.3321893233.all_goals_core._main._lambda_2
7388ms    98.7%   all_goals_core
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 (Mar 15 2019 at 15:02):

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 (Mar 15 2019 at 15:06):

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 (Mar 15 2019 at 15:28):

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

Last updated: Aug 19 2019 at 18:00 UTC