Discussion: ComputableReal

Stream: general

[Alex Meiburg]

Making a thread for ComputableReal, in case anyone has questions about it. :)

[Rob Lewis]

This is very cool -- could you write a little about the representation you're using?

[Yaël Dillies]

Nice! It looks like you're using docs#NonemptyInterval, which makes me feel a bit proud :blush:

[Yaël Dillies]

Your instance here duplicates @Eric Wieser's instance from #20636. Would you be interested in getting the PR over the line?

[Alex Meiburg]

@Rob Lewis I just finished writing a Readme that I hope helps: https://github.com/Timeroot/ComputableReal/blob/main/README.md

[Alex Meiburg]

@Yaël Dillies : Sure! I tried proving that multiplication was associative, though, and it kind of blew up. Just a whole lot of different cases and not an obvious easy way to clean things up. So I don't know how easy getting the Semiring instance on intervals would be.

[Vlad Tsyrklevich]

As a good practice, I’d also include a license, probably Apache to match the ecosystem

[Kevin Buzzard]

I know nothing about this sort of thing, but would it be possible to make decide work? I can't use theorems which are only proved using native_decide because it's easy to prove False using it.

[Alex Meiburg]

Yes, it should be! Right now there are a couple of 'bad decidability instances' in Mathlib that don't play well decide (but native_decide isn't so picky), but that should be a pretty quick PR.

[Alex Meiburg]

No promises that it will be particularly efficient. The code was written to be efficient when compiled, but I have no idea how decide will do with it.

[Kevin Buzzard]

But most real numbers of interest can be computed.

In my line of work, many natural numbers of interest cannot be computed! They are "the number of solutions to this polynomial equation" and there are no known (provably correct) algorithms for answering questions such as this. For example there are certainly only finitely many rational solutions to $x^4+y^4=17$ by a deep theorem of Faltings, but figuring out how many there are is a very different game.

[Bolton Bailey]

I am trying to test this out (this was a goal I needed for Carleson a few weeks ago), but I am getting the following error message

import ComputableReal.ComputableRSeq
import ComputableReal.ComputableReal
import ComputableReal.IsComputable
import ComputableReal.SpecialFunctions

example : Real.sqrt 2 > 4 / 3 := by --true inequality, good
  native_decide -- works

lemma geom_estimate_constant_le_two :
    1 ≤ (2 : ℝ) * (4 * (1 - 2 ^ (-1 / 4 : ℝ))) := by
  native_decide -- failed to compile definition, consider marking it as 'noncomputable' because it depends on 'Real.instDivInvMonoid', and it does not have executable code

[Alex Meiburg]

I would say that all those real numbers can be computed, but we simply don't know how! Certainly for any natural number, there's an algorithm that computes it. That might feel pedantic, but I do think that's the standard meaning of "computable number" as it shows up in complexity theory.

I understand your intent, though... like there's some function that takes certain integer data (like coefficients of some polynomial equation) and gives you a natural number (the number of solutions) that are not known to be computable. And in some cases, of course, is known not to be computable.

[Bolton Bailey]

Ah, perhaps negative powers simply aren't supported yet? (A more careful reading of the README on my part tells me this is the case)

[Alex Meiburg]

Bolton: rational powers aren't supported at the moment. :/ But negative integers are

[Kevin Buzzard]

If that's what computed means then is it possible to give any example of a real number which can't be computed in that sense?

[Alex Meiburg]

Sure, like Chaitin's constant. The existence of an algorithm that computes a convergent sequence of bounds on that constant, would necessarily also give you an algorithm to solve the halting problem.

[Bolton Bailey]

Regardless of the feature set, this is awesome! I have wanted something like this since at least 2021. Are you accepting contributions?

[Alex Meiburg]

Gladly! :grinning:

[Eric Wieser]

Is replacing mathlib's Real with your version a viable option in the longer term?

[Eric Wieser]

(if I understand correctly, instead of storing one cauchy sequence, you store two where one is pointwise less than or equal the other?)

[Edward van de Meent]

Kevin Buzzard said:

If that's what computed means then is it possible to give any example of a real number which can't be computed in that sense?

yes. since there are countably many possible computations, and uncountably many reals, there must be uncountably many that aren't represented by a computation.

[Mario Carneiro]

When I first saw this I found it a bit unbelievable, since it seems to claim that real number equality is decidable by reduction to computable reals, but equality of computable reals is also undecidable. Turns out the trick is that it's not actually a decision procedure, there is a way to trick the compiler into compiling a nonterminating function for the decision procedure and thereby getting a lean-computable Decidable instance which isn't a decision procedure

[Mario Carneiro]

But I agree with Eric, we should make this the actual definition of reals, because it will allow us to make lots of things computable which is otherwise a bit of a pain to work around and a new keyword to learn about

[Alex Meiburg]

Eric Wieser said:

(if I understand correctly, instead of storing one cauchy sequence, you store two where one is pointwise less than or equal the other?)

Almost - just being pointwise less isn't enough (indeed I made that mistake at first), the lower bounds need to be less than the point they converge to. For instance, a sequence of lower bounds on pi couldn't include 3.2, even if the upper bound at that point in the sequence was 3.5.

[Mario Carneiro]

It should be possible to require that the sequences are monotonic

[Mario Carneiro]

Do your sequences have a modulus of convergence?

[Mario Carneiro]

I'm inclined to think they should

[Eric Wieser]

Mario Carneiro said:

When I first saw this I found it a bit unbelievable, since it seems to claim that real number equality is decidable by reduction to computable reals, but equality of computable reals is also undecidable. Turns out the trick is that it's not actually a decision procedure, there is a way to trick the compiler into compiling a nonterminating function for the decision procedure and thereby getting a lean-computable Decidable instance which isn't a decision procedure

I guess this is an argument for the semidecision procedures discussed elsewhere?

[Mario Carneiro]

I think that wouldn't help here

[Mario Carneiro]

you can't ever return none from functions that observe the cauchy sequences because it will usually have to depend on the details of the sequence, so it won't be stable under the quotient

[Alex Meiburg]

Yeah, you can require monotonicity, and you can also make a sort of wrapper to turn any non-monotonic form into a monotonic form by just keeping the best bound found so far. The reason I didn't do this in this package was

• If you require monotonicity, then you also to prove that the sequences are monotonic at each step, and that's a considerable extra burden
• If you don't require it, but do the wrapping, then you get a nonnegligble extra computational overhead ... unless you're careful+smart about how you store that running bound.

[Alex Meiburg]

Mario Carneiro said:

But I agree with Eric, we should make this the actual definition of reals, because it will allow us to make lots of things computable which is otherwise a bit of a pain to work around and a new keyword to learn about

I'm surprised to hear you both say this! Personally I don't think it's very viable or useful to replace Mathlib's Real with this in the long term. It's possible in principle - one can classically prove that "all reals are computable", in the same way that "all propositions are Decidable" - but I don't know what it would get you, really. I also think that, architecturally, _first_ defining reals (so that you have the nice existence and topology and theory of them) and _then_ talking about which ones can be computed, works out better.

[Mario Carneiro]

If you have a fixed rate of convergence, I think you can get away with just having one sequence instead of two

[Mario Carneiro]

and you get monotonicity for free

[Edward van de Meent]

Kevin Buzzard said:

If that's what computed means then is it possible to give any example of a real number which can't be computed in that sense?

looking into it a bit, Chaitin's constant is one among several "concrete" examples of noncomputable reals

[Alex Meiburg]

Yeah, I definitely don't want to have a requirement on rate of convergence, since there are some computable functions that don't "naturally" have fast convergence -- so then you need something where you accelerate the sequence (I think that's the technical term for it) to get the desired rate

[Mario Carneiro]

that is, it's a sequence of integers a : Nat -> Int such that x * 2^n is within 1 of a n

[Mario Carneiro]

Alex Meiburg said:

Yeah, I definitely don't want to have a requirement on rate of convergence, since there are some computable functions that don't "naturally" have fast convergence -- so then you need something where you accelerate the sequence (I think that's the technical term for it) to get the desired rate

This is another trivial wrapper operation, like the monotonicity one you mentioned

[Mario Carneiro]

Just to be clear: this definition does not give computable reals, it gives all reals. The computable reals are the ones where you can write down the definition of the real number without using noncomputable (assuming you don't cheat)

[Alex Meiburg]

Yes, maybe I should make that more clear in the Readme, that it doesn't let you actually prove computability or anything.

[Mario Carneiro]

So it would work out fine as a replacement for mathlib reals, at least from the perspective of provable theorems

[Mario Carneiro]

and I do like the idea that you can do computations on reals and run #eval and get an answer that isn't useless

[Alex Meiburg]

I think there's a natural sort of tension between "definition of reals that lets you prove theorems easily" and "definition of reals that lets you compute things quickly". Of course with enough theorems, you can reproduce the current Real API 100%, so the first point would become moot. So I'm still skeptical

[Kevin Buzzard]

Edward van de Meent said:

Kevin Buzzard said:

If that's what computed means then is it possible to give any example of a real number which can't be computed in that sense?

yes. since there are countably many possible computations, and uncountably many reals, there must be uncountably many that aren't represented by a computation.

I'm well aware that there must exist noncomputable reals, I'm just saying that when you're talking about a specific real number rather than a generic one (I'm interpreting a "real number of interest" as "not an arbitrary real number") then somehow it must be computable by your definition, because the very fact that you're talking about it means that it's somehow been uniquely expressed in your system.

[Alex Meiburg]

I think that's the notion of "definable real", which is a countable strict superset of the computable reals.

[Edward van de Meent]

that's not true. see above mentioned "Chaitin's constant"

[Edward van de Meent]

also, numbers like "the real number whose binary representation encodes the solution to the halting problem" and variants

[Edward van de Meent]

those are definitely definable, but not computable

[Michael Stoll]

An implementation of computable reals in Common Lisp can be found here (this is based on code I wrote a long time ago...). From the README:

Computable real numbers x are interpreted as (potentially) infinite fractions in base 2 that are specified through a rule for computation of an integer a with |(2^k)*x - a| <= 1 for any k>=0.

This is quite close to what Mario wrote above. (This is, of course, a bit off-topic here, but may provide some inspiration.)

[Mario Carneiro]

I'm not sure I would say it is off topic! Formalization is at least 70% imitation of things that have been done before in a different way

[Mario Carneiro]

I think there is also a computable reals library in haskell somewhere, it's probably worth reviewing existing libraries if we want to find a good representation

[Rob Lewis]

There's a long line of research into computable representations of the reals and no obvious "best" option, it really depends on your applications. My memory is that Cauchy sequences with a fixed modulus of convergence are a compromise that allow for a little computation in practice, without horribly complicating your proofs.

[Rob Lewis]

I'm no expert here, but there are pointers to a bunch of different approaches in the proceedings of a CIRM meeting I went to ages ago.

[Mario Carneiro]

One nice thing that the decidability trick adds to the equation is the ability to actually get digit representations out of the quotient, so this would eliminate the need for our unsafe Repr Real instance

[Alex Meiburg]

Kevin Buzzard said:

I know nothing about this sort of thing, but would it be possible to make decide work? I can't use theorems which are only proved using native_decide because it's easy to prove False using it.

Hmm, maybe not so easy. I thought it was a bad decidability instance -- it wasn't, my mistake. I think the actual issue is that ComputableRSeq.sign is (by necessity) a partial def. It seems that decide and decide +kernel won't look inside those. The docs say,

Decidable instances defined by well-founded recursion might not work because evaluating them requires reducing proofs.

which hadn't clicked for me before partial isn't well-founded recursion. Well, it's not-well-founded recursion, which I suppose makes it strictly worse!

[Mario Carneiro]

You are right, this trick will not work in the kernel. However, you can still do proofs using it using a slightly different interface, which takes a fuel amount as input and produces a maybe-proof, and if the fuel is high enough then you get your proof (and a tactic can hide the details of that)

[Mario Carneiro]

that was always how I wanted the approx tactic to work

[Niels Voss]

Discussion on encoding possibly-non-terminating computations in an optimal way is happening here: #Is there code for X? > Divergence monad

[Bas Spitters]

Eric Wieser said:

Is replacing mathlib's Real with your version a viable option in the longer term?

This is the current design in Coq, and generally the philosophy of Bishop's mathematics. By just adding excluded middle, one obtains the classical results.

[Bas Spitters]

@Alex Meiburg I'm curious how your timings compare to the ones we obtained in MathClasses(https://math-classes.github.io/). Specifically: https://lmcs.episciences.org/958
Is this was 12 years ago, you should get a good speed-up due to Moore's law.

[Jihoon Hyun]

I'm happy to see someone's working on this! (I was thinking about this topic a lot, but never started it.)
The fastest non-lean library I know which computes reals is iRRAM. Maybe one can compare or get an idea from this library?
Also some reference: Computable Analysis: An Introduction

[Bas Spitters]

We tried to explain some of the issues with real computation and non-computable numbers here

[Bas Spitters]

@Jihoon Hyun the main observation in iRRAM is that one does forward computation. One tries single precision interval arithmetic. If this in not good enough, one starts over with double precision. The usual Cauchy representation computes backwards.

[Bas Spitters]

@Jihoon Hyun Competition between the various systems:
https://www.cs.ru.nl/~freek/pubs/comp_rep.pdf

[Bas Spitters]

@Alex Meiburg You should be able to get a substantial speed boost by using dyadic rationals. (at least that was the case for us).

Your approach using an IsComputable type class seems similar that computing with classical reals in Coq.

[Alex Meiburg]

Mario Carneiro said:

You are right, this trick will not work in the kernel. However, you can still do proofs using it using a slightly different interface, which takes a fuel amount as input and produces a maybe-proof, and if the fuel is high enough then you get your proof (and a tactic can hide the details of that)

I thought about this some more. Yeah, doing it with a tactic is always one way to go about it, but I am quite keen on making this work with decide in particular since then it's part of this nice system. I think it can be made to work with fuel in this kind of way -- something like (getting of pseudocode-y here)

--Partial decidable statements: you may or may not get a decidable witness out
class PDecidable (p) := Option (Decicdable p)

instance PDecidableLE : PDecidableLE ComputableReal :=
   checkSignWithFuel 100000 --gives `some Decidable` if we successfully proved/disproved it, otherwise `none`

noncomputable instance Decidable_of_PDecidable {p : Prop} (pd : PDecidable p) : Decidable p :=
   @decidable_of_iff _ (if h : Option.isSome pd then p else p) (by rw [dite_eq_ite, ite_self])
    (instDecidableDite (dT := fun h ↦ Option.get _ h) (dE := fun _ ↦ Classical.propDecidable _))

So then when decide runs, it will first evaluate the PDecidable with the given fuel, if it works and gets some Decidable then it uses that, otherwise it falls back to Classical.propDecidable and then gives up and fails. I tried this out on some toy examples, this seems to work with decide just fine!

[Alex Meiburg]

Of course it doesn't work with native_decide, since that requires everything to be computable at compile-time, and Classical.propDecidable isn't. And because there is no complete decision procedure for reals, you can't have an instance that works both with decide and native_decide.

I think making each of these decidability instances a separate scoped instance is a good solution then. (I also, in general, think that it should be scoped anyway, because having a not-necessarily-terminating instance "running wild" does make me a bit uncomfortable.)

[Alex Meiburg]

So then it would be something like

theorem : .... := by
   decide

and

theorem : .... := by
   native_decide

and they should both work.

[Niels Voss]

Just a thought, but could you replace the Classical.propDecidable in your instance with a partial def? From the perspective of the kernel, partial def is pretty similar to something produced by choice, but it lets the compiler keep going.

[Niels Voss]

You would use Classical.propDecidable to prove that such a decidability instance actually exists (since partial def needs a Nonempty instance) but you don't need to mark anything noncomputable

[Alex Meiburg]

That's true, I could do that! Or in other words, just make the native_decide instance be the fallback.

[Alex Meiburg]

Then decide will use up the fuel, and if it hasn't determined an answer, hit the partial def and quit. But native_decide will be free to use the partial def and keep going arbitrarily long (at the slight expense of redoing the first fuel many steps of the process, maybe)

[Mario Carneiro]

it won't quit at the partial def, it will loop

[Alex Meiburg]

Why would it loop? decide sees an opaque def and gets stuck

[Mario Carneiro]

right, sorry for the confusion. So here's a toy version of your program then:

partial def findTrue (f : Nat → Bool) (i : Nat) : Decidable (∃ i, f i) :=
  if h : f i then .isTrue ⟨_, h⟩ else findTrue f (i + 1)

def findTrueWithFuel (f : Nat → Bool) (i : Nat) : Nat → Decidable (∃ i, f i)
  | 0 => findTrue f i
  | fuel + 1 => if h : f i then .isTrue ⟨_, h⟩ else findTrueWithFuel f (i + 1) fuel

[Mario Carneiro]

and the latter function works in both the kernel and the compiler if you give it enough fuel

[Alex Meiburg]

Yeah, I think that's an accurate picture

[Mario Carneiro]

and of course a good compiler should notice that the function can only ever return isTrue and compile it to just return true :upside_down:

[Alex Meiburg]

Haha, in the toy example, yes. :) But in my case it can return true or false. (If the lower and upper bounds are both equal to 0 at some point, then we know the true real value is zero, and so can we prove equality.)

[Niels Voss]

Wait, if the compiler would make it return true instead of entering an infinite loop, wouldn't that make it easy to prove False with native_decide by taking a sequence of all falses and showing it must be true at some point?

[Mario Carneiro]

infinite loops are UB in C, which allows for some very powerful (or unsound, depending on who you ask) optimizations

[Mario Carneiro]

and lean compiles to C, so...

[Alex Meiburg]

Good news, it does not seem to let you use native_decide to prove False .... at the moment. (I agree, it seems tenuous, like a slight change in compiler could mess it up.) Bad news: you can write fine code with no sorries and then running #eval throws an unreachable code error.

import Mathlib

partial def isZeroOne (_ : ℕ) : { b : Bool // b = (1 == 0) } :=
  if h : 1 = 0 then ⟨true, h⟩ else isZeroOne 0

instance decZeroOne : Decidable (1 = 0) :=
  decidable_of_iff (isZeroOne 0) ⟨(by simp [(isZeroOne 0).2] at ·), by simp⟩

#eval decZeroOne --unreachable code

theorem zero_is_one : 1 = 0 := by
  -- native_decide --hangs

(Tried this in Lean4 web, live)

[Mario Carneiro]

GCC gets close, but does not seem to perform the optimization I was hoping to see.

lean:

partial def foo (b : Bool) := if b then true else foo b

C:

uint8_t l_foo(uint8_t x_1) {
_start:
{
if (x_1 == 0)
{
goto _start;
}
else
{
uint8_t x_3;
x_3 = 1;
return x_3;
}
}
}

GCC -O2:

l_foo:
        cbz     r0, .L3
        movs    r0, #1
        bx      lr
.L3:
        b       .L3

The last two lines are the infinite loop. It got optimized to if b then true else loop() while I was hoping it would just be true, which should be legal per the C spec

[Niels Voss]

If you use #eval is the compiler being invoked, or does it use some other sort of code runner? I thought I read somewhere that there were a couple different ways to run Lean code, but maybe I'm misremembering.

As for the C code, maybe try turning up to max compiler optimizations to see if it performs the loop optimization then?

[Mario Carneiro]

#eval uses the lean compiler pipeline as far as generating optimized IR, then runs the interpreter

[Mario Carneiro]

above I'm showing the result of compilation to C, which you won't be running if you use #eval

[Mario Carneiro]

the lean compiler is actually quite conservative around partial defs because it has been hit in the past by bugs exactly like the ones we're talking about

[Niels Voss]

Alex Meiburg said:

Good news, it does not seem to let you use native_decide to prove False .... at the moment. (I agree, it seems tenuous, like a slight change in compiler could mess it up.) Bad news: you can write fine code with no sorries and then running #eval throws an unreachable code error.

import Mathlib

partial def isZeroOne (_ : ℕ) : { b : Bool // b = (1 == 0) } :=
  if h : 1 = 0 then ⟨true, h⟩ else isZeroOne 0

instance decZeroOne : Decidable (1 = 0) :=
  decidable_of_iff (isZeroOne 0) ⟨(by simp [(isZeroOne 0).2] at ·), by simp⟩

#eval decZeroOne --unreachable code

theorem zero_is_one : 1 = 0 := by
  -- native_decide --hangs

(Tried this in Lean4 web, live)

From the docstring of lcUnreachable:

Executing this expression to actually synthesize a value of type α causes immediate undefined behavior, and the compiler does take advantage of this to optimize the code assuming that it is not called. If it is not optimized out, it is likely to appear as a print message saying "unreachable code", but this behavior is not guaranteed or stable in any way.

[Niels Voss]

So unfortunately, it is indeed UB

[Mario Carneiro]

I wrote that line :)

[Mario Carneiro]

Mario Carneiro said:

the lean compiler is actually quite conservative around partial defs because it has been hit in the past by bugs exactly like the ones we're talking about

but I suspect this is not good enough because it generates infinite loops without any operations in them in generated C code, so the C compiler that runs on this will not necessarily be similarly conservative

[Mario Carneiro]

but it seems Alex found a bug in the lean compiler anyway

[Alex Meiburg]

I'm more than slightly proud of this. :)

[Niels Voss]

I am just now realizing how even the slightest mistake in an @implemented_by clause (which I think you or someone else called a "compiler axiom") can lead to UB. I always knew it could cause the proof semantics to mismatch the program semantics, but I didn't realize it could lead to UB

[Niels Voss]

I probably should have known though because unwrapping a Part can lead to UB, as discussed in the other thread linked above

[Mario Carneiro]

here's an easy example of UB from any kind of divergence between compiler semantics and proof:

def liarImp := true
@[implemented_by liarImp]
def liar := false

#eval #[0][if liar then 1000 else 0] -- 727562280

[Niels Voss]

Is this related to Array being defined outside Lean? Does it still cause UB for List?

[Mario Carneiro]

it manifests differently, but both are UB

[Mario Carneiro]

Array makes it easy to do arbitrary buffer overruns

[Mario Carneiro]

with List you unwrap an option that isn't an option resulting in a pointer dereference of something that isn't a pointer

[Mario Carneiro]

which you can also weaponize with some ingenuity

[Alex Meiburg]

Those both require implemented_by or unsafe calls, right?

[Mario Carneiro]

or compiler bugs

[Alex Meiburg]

Right haha

[Mario Carneiro]

anything leading to some decidable property which evaluates to true and is provably false

[Mario Carneiro]

not having any such thing is a really strong constraint on lean code, it means for example that you can't convert a Float to a UInt64 without masking out nans

[Joseph Myers]

Infinite loops are UB in C++. The forward progress guarantees in C are much weaker than in C++ (only for iteration statements whose condition is not a constant expression, not for loops with goto or with constant true controlling the infinite loop).

[Mario Carneiro]

You're right. If I compile the same code as C++ instead of C the loop disappears:

_Z5l_fooh:
        movs    r0, #1
        bx      lr

this is just return true

[Mauricio Collares]

(deleted)

[Alex Meiburg]

I just learned that Mathlib already has a Estimator for doing something very similar, while caching some data.

Perhaps some day it would be nice to rewrite this to use that, which would definitely allow for some significant runtime improvements.

[Kim Morrison]

Estimator is very specialized for its use case inside rw_search, and it's unclear that it's a good design (mine!) for general purpose use.

[Alex Meiburg]

Well, it is pretty much exactly (up to any Lean-specific idioms I do or don't know) the design I had imagined for a while: a value, and (in particular) a Thunk for improving it, with kind of a black box as to how much it improves.

But, don't worry, I won't be attempting to use it any time soon. :) I'm working on getting log supported, and that's giving me quite a time already.

[Alex Meiburg]

And always a zillion other things..