Computable real.sqrt

Stream: maths

[Patrick Stevens]

If anyone is interested in a mini project to fill in a long-standing TODO in data.real.sqrt, I've started https://github.com/leanprover-community/mathlib/pull/13634 which adds some theorems around sqrt_aux, in preparation for showing that Heron's method defines a sequence which is Cauchy for sufficiently nicely-expressed input reals, and then for showing that it defines a real number which is equal to the square root of the input (I've already got that it's at least the square root). Feel free to pitch in!

I'm particularly sad about the amount of grubbing around I had to do to show the most basic things about rational numbers.

[Kevin Buzzard]

A lot of the grubbing around is stuff we have already. Do you know about library_search? I left some comments.

[Patrick Stevens]

Thanks - I used library_search extensively to get that far, and I have no idea why I didn't find those, but the ones you pointed out certainly work when I run library_search now! I thought I was being careful only to write a proof if I really couldn't find it, and I certainly think that library_search is the first line I write when I try and prove something; I'll put my bizarre failure down to my new and shiny covid infection.

[Patrick Stevens]

Sigh, I really should learn to run little experiments before trying to prove something - I just realised I spent an awful lot of effort trying to prove something that's patently false

[Patrick Stevens]

But I did get an epsilon down quite small - it's just not an epsilon that's any use to anyone

[Patrick Stevens]

Does anyone know whether Heron's formula actually converges at 0?

def sqrt_aux (f : cau_seq ℚ abs) : ℕ → ℚ
| 0       := max 1 (f 0)
| (n + 1) := let s := sqrt_aux n in (s + (max 0 (f (n+1))) / s) / 2

I've shown that it converges and is equal (as a real number) to sqrt f if f > 0 or if f < 0, and I've shown that if f >= 0 then the sequence infinitely often gets arbitrarily close to the square root, but for the life of me I can't show that it converges when f = 0 (as a real number). Does anyone know if it's even true?

[Patrick Stevens]

(a number of useful facts, like sqrt_aux f i ^ 2 >= f i for all i, are present in https://github.com/leanprover-community/mathlib/pull/13634/files)

[Patrick Stevens]

Defining f[n_] := If[OddQ[n], 0, 1/2^(n/2)] as a representation of 0, and using an initial guess of 0.01 to simulate a difficult case, Mathematica gives {0.005, 50.0025, 25.0013, 12.5056, 6.25281, 3.1364, 1.5682, 0.804028, 0.402014, 0.239874, 0.119937, 0.125107, 0.0625534, 0.0937233, …, which certainly doesn't seem to be converging in any nice-and-easy-to-prove way

[Reid Barton]

Looks to me like it does not converge. If f n is 0 for many successive n then the approximations can get arbitrarily small, and then if the next f n is size ~epsilon, the next approximation could be unboundedly large

[Patrick Stevens]

I have increasingly been coming to that conclusion :(

[Reid Barton]

What I think would work is to first do this for a fixed rational number, and then diagonalize.

[Reid Barton]

But anyways, cau_seq is not a reasonable or useful notion of computability for real numbers

[Patrick Stevens]

All I really wanted to do was remove the word noncomputable from the current real.sqrt

[Patrick Stevens]

Do you mean cau_seq is not useful because it's just an existential quantifier which is a Prop rather than a proper sigma type with data in?

[Reid Barton]

noncomputable theory

[Reid Barton]

Right, there is no information about the rate of convergence

[Jireh Loreaux]

@Patrick Stevens If I give you a Cauchy sequence and tell you the billionth term is 5, you can tell me exactly zero information about the real number corresponding to that Cauchy sequence.

[Patrick Stevens]

Well, only because our "Cauchy" doesn't tell you the values inside the existentials, but yes

[Jireh Loreaux]

Right, you would need something like Bishop's constructive reals

[Patrick Stevens]

Ah, given that it's impossible to compute an $A$ out of an $A / \sym$, this was never going work in the first place (as I'm sure was obvious to nearly everyone) - I can't get a cau_seq out of an $\mathbb{R}$ at all

[Patrick Stevens]

(I'm still not used to Prop :P I still expect things to have data attached to them)

[Eric Wieser]

It's possible to do that extraction in meta land with docs#quot.unquot

[Reid Barton]

This wasn't going to be a problem, you can use quotient.lift as long as the resulting real doesn't depend on the choice of representative (which must be the case if it was going to be the square root)

[Kyle Miller]

Given some semi-recent changes to Lean 3, you can define an #eval-friendly typeclass that gives how to compute particular reals using rational approximations:

class real.computable (x : ℝ) :=
(approx (ε : ℚ) (pos : 0 < ε) : ℚ)
(is_close (ε : ℚ) (pos : 0 < ε) : |(approx ε pos - x : ℝ)| < ε)

I defined some instances for bit0, bit1, some coercions, and addition. Here's an approximation of 2+4 that's within 1/10:

#eval (2 + 4 : ℝ).approx (1/10) (by linarith)
-- 6

code

import data.real.basic
import tactic

class real.computable (x : ℝ) :=
(approx (ε : ℚ) (pos : 0 < ε) : ℚ)
(is_close (ε : ℚ) (pos : 0 < ε) : |(approx ε pos - x : ℝ)| < ε)

alias real.computable.approx ← real.approx

lemma real.approx_is_close (x : ℝ) [x.computable] (ε : ℚ) (pos : 0 < ε) :
  |(x.approx ε pos - x : ℝ)| < ε :=
real.computable.is_close ε pos

instance real.computable.nat (n : ℕ) : real.computable n :=
{ approx := λ ε pos, n,
  is_close := by simp }

instance real.computable.has_zero : real.computable 0 :=
{ approx := λ ε pos, 0,
  is_close := by simp }

instance real.computable.has_one : real.computable 1 :=
{ approx := λ ε pos, 1,
  is_close := by simp }

instance real.computable.int (n : ℤ) : real.computable n :=
{ approx := λ ε pos, n,
  is_close := by simp }

instance real.computable.add (x y : ℝ) [x.computable] [y.computable] : (x + y).computable :=
{ approx := λ ε pos,
    x.approx (ε / 2) (half_pos pos) + y.approx (ε / 2) (half_pos pos),
  is_close := λ ε pos, begin
    have := half_pos pos,
    have hx := x.approx_is_close (ε / 2) ‹_›,
    have hy := y.approx_is_close (ε / 2) ‹_›,
    calc |↑(x.approx (ε / 2) ‹_› + y.approx (ε / 2) ‹_›) - (x + y)|
        = |(↑(x.approx (ε / 2) ‹_›) - x) + (↑(y.approx (ε / 2) ‹_›) - y)| : _
    ... ≤ |↑(x.approx (ε / 2) ‹_›) - x| + |↑(y.approx (ε / 2) ‹_›) - y| : abs_add _ _
    ... < ↑(ε / 2) + ↑(ε / 2) : add_lt_add hx hy
    ... = ↑ε : by simp,
    { congr,
      norm_cast,
      simp [← sub_sub, sub_eq_add_neg, add_assoc, add_left_comm, add_comm], },
  end }

instance real.computable.bit0 (x : ℝ) [x.computable] : (bit0 x).computable :=
{ approx := λ ε pos, let q := x.approx (ε / 2) (half_pos pos) in 2*q,
  is_close := λ ε pos, begin
    simp [two_mul],
    let := real.computable.add x x,
    -- there is a bug in dot notation (patch in progress)
    let := this.is_close,
    convert this ε pos,
    simp [computable.approx],
  end }

instance real.computable.bit1 (x : ℝ) [x.computable] : (bit1 x).computable :=
{ approx := λ ε pos, let q := x.approx (ε / 2) (half_pos pos) in 2*q + 1,
  is_close := λ ε pos, begin
    have : ∀ y, y + 1 - bit1 x = y - 2*x,
    { simp [bit1, bit0, add_mul], },
    simp [two_mul, this],
    let := real.computable.add x x,
    -- there is a bug in dot notation (patch in progress)
    let := this.is_close,
    convert this ε pos,
    simp [computable.approx],
  end }

#eval (2 + 4 : ℝ).approx (1/10) (by linarith)
-- 6

[Kyle Miller]

These sorts of instances are a place you could put an algorithm for square roots of reals.

[Kyle Miller]

By the way, that real.computable class is just proof-of-concept. It seemed plausible as a definition, but there could be better ones.

[Mario Carneiro]

This is interesting, because I recently investigated this question with regard to a proof in intuitionistic logic of the existence of the real square root function using Heron's formula in metamath, and the proof does go through: The sequence $x_{n+1}=(x_n+a/x_n)/2$ converges to $\sqrt a$ for all $a\ge 0$.

But @Reid Barton is right that this proof does not work when you use a Cauchy sequence $a_n$ instead of literally using $a$ there. Concretely, we can define a sequence $a_n$ such that $a_0=1$ and for other $n$, if $n=k^2$ is a perfect square then $a_n=8(1-2/4^k)/4^k$, and otherwise $a_n=0$. Then sqrt_aux a has the pattern $1,1,2^{-1},2^{-2},1,2^{-1},2^{-2},2^{-3},2^{-4},1,\dots$: it gets arbitrarily close to $0$ and is also equal to $1$ infinitely often, so it does not converge. But $0\le a_n\le 8/4^k$ for $n>0$ so it converges to 0.

[Patrick Stevens]

When it's less late at night I'll play with sqrt_aux x e n having max (e ^ 2) (f (n + 1)) instead of the current e = 0, and see if I can make that converge, and then hope that sqrt x n := sqrt_aux x (1 / n) n does the job

[Mario Carneiro]

Do we have the ability to make computable cauchy sequences of reals? That would be the easiest way to make the standard proof work