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


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
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)) _,
norm_num
end


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

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

Johan Commelin (Mar 15 2019 at 04:53):

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

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/

Kevin Buzzard (Mar 15 2019 at 09:31):

Comments welcome as ever.

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
It doesn't require login

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
# rule: 2+x<=y^2 |- sqrt_two_add_series(x,n+1) <= sqrt_two_add_series(y,n)

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:18):

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

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

Kenny Lau (Mar 15 2019 at 11:39):

I already optimized this

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