Numerical approximations of real valued functions

Stream: Is there code for X?

[Niels Voss]

Is there a library which for some functions f : Real -> Real, lets you approximate f with a computable g : Float -> Float (perhaps by recursively looking at the definition of f and replacing Real.sin with Float.sin, Real.pi with a floating point pi, etc.)? I had thought that maybe SciLean was able to do this, but after I did a bit of looking into it it didn't seem like this was one of SciLean's capabilities.

[Aaron Liu]

see related #mathlib4 > Numeric estimates for logarithms

[Niels Voss]

ComputableReal looks promising. I don't actually need proofs, and it's fine if it's theoretically possible for numerical inaccuracies to build up and make the results unusable, as long as it doesn't happen all the time

[Niels Voss]

The reason I'm interested in this is that I'm trying to see if it is possible to make plausible work on expressions with real numbers by evaluating using floats and checking if their difference exceeds some predefined epsilon

[Jireh Loreaux]

I think the answer to your latter query is no, because you can't reason about Float, it's opaque.

[Niels Voss]

Oh I just looked at the API for plausible; I didn't realize that it required you to actually provide a proof that the proposition is false when you provide a counterexample

[Niels Voss]

Plausible uses native_decide when providing counterexample proofs, right?

[Niels Voss]

import Plausible

opaque n : Nat := 37
#test n = 37 -- Unable to find a counter-example
#test n = 22 -- Found a counter-example!

[Niels Voss]

So I guess this means that ComputableReal might work

[Antoine Chambert-Loir]

Jireh Loreaux said:

I think the answer to your latter query is no, because you can't reason about Float, it's opaque.

Does that mean that we couldn't prove in Lean algorithms involving docs#Float ?

[Jireh Loreaux]

If you wanted to reason about Float, you could, but you would need to add axioms about it's behavior. This is rather dangerous because the IEEE spec has many bad properties, and I think many implementations don't conform to the spec on all points. I'm speaking from memory here.

[Alex Meiburg]

"someone should" write down the axioms. The most I'm aware of is https://github.com/bustercopley/lean-float .

But yes, if you're talking about "floats" and "proofs" in the same sentence, and it's for any reason besides "I need 100% reproducible scientific computing", you've probably done something wrong.

[Kenny Lau]

(from brief reading, that lean-float doesn't seem to be exactly the same as, well, the IEEE float)

[Dominic Steinitz]

I always find it helpful with Float to imagine a world in which e.g. the exponent is 2 bits and the mantissa is 3 bits. Now you can really see the issues with using Float.

I didn't get very far with this but I did try using a type class with Float and Rational as members. One of the problems I ran into was what is $\sqrt{2}$. I decided on some approximation for the Rationals and my type class had a method (sorry for using Haskell terminology) for $\sqrt{2}$. But beware: https://idontgetoutmuch.wordpress.com/2015/11/12/floating-point-a-faustian-bargain/

Also there is literature on this subject, e.g.: https://icfp18.sigplan.org/details/npfl-2018-papers/5/Exact-Real-Arithmetic-for-Geometric-Operations

[Jz Pan]

For pure software IEEE754 float implementation in Lean discussion, maybe you can see #Is there code for X? > Proving nan != nan @ 💬

[Niels Voss]

But I don't actually need proofs; it's fine if there are some false positives

[Niels Voss]

I guess what I'm looking for is something like this:

import Mathlib.Data.Real.Basic
import Mathlib.Data.Complex.Trigonometric
import Mathlib.Analysis.SpecialFunctions.Trigonometric.Basic
import Std
import Qq

open Qq Lean

def Float.pi : Float := 3.1415926535897931

partial def convert (e : Q(ℝ)) : MetaM (Q(Float)) :=
  match e with
  | ~q(OfNat.ofNat 0) => return q(0 : Float)
  | ~q(OfNat.ofNat 1) => return q(1 : Float)
  | ~q(@OfNat.ofNat ℝ $x (@instOfNatAtLeastTwo ℝ $x Real.instNatCast $p)) =>
    return q(Float.ofNat $x)
  | ~q(@OfScientific.ofScientific ℝ NNRatCast.toOfScientific $mantissa $negative $exponent) =>
    return q(Float.ofScientific $mantissa $negative $exponent)
  | ~q(-$x) => return q(-$(← convert x))
  | ~q($x + $y) => return q($(← convert x) + $(← convert y))
  | ~q($x - $y) => return q($(← convert x) - $(← convert y))
  | ~q($x * $y) => return q($(← convert x) * $(← convert y))
  | ~q($x / $y) => return q($(← convert x) / $(← convert y))
  | ~q(Real.sqrt $x) => return q(Float.sqrt $(← convert x))
  | ~q(Real.sin $x) => return q(Float.sin $(← convert x))
  | ~q(Real.cos $x) => return q(Float.cos $(← convert x))
  | ~q(Real.tan $x) => return q(Float.tan $(← convert x))
  | ~q(Real.sinh $x) => return q(Float.sinh $(← convert x))
  | ~q(Real.cosh $x) => return q(Float.cosh $(← convert x))
  | ~q(Real.tanh $x) => return q(Float.tanh $(← convert x))
  | ~q(Real.pi) => return q(Float.pi)
  | _ => throwError f!"could not convert"

elab "asdf" : term => convert q(Real.sqrt (-4 * Real.cos (-Real.pi)) : ℝ)

#eval asdf -- 2.000000

[Niels Voss]

but more robust, and maybe using typeclasses instead of pattern matching

[Richard Copley]

Kenny Lau said:

(from brief reading, that lean-float doesn't seem to be exactly the same as, well, the IEEE float)

Indeed, nowhere near. Just enough for some results about addition and cut-subtraction of positive floating-point numbers, in certain rather idiosyncratic classes of floating-point arithmetics, assuming overflow does not occur [edit: in the 'exponent' field]. And there is no Float type there; the results need some interpretation.

[Tyler Josephson ⚛️]

Is this approach of interest?
Our typeclass is called RealLike, and contains the instances we need for our LeanLJ project. You can go ahead and add additional instances for the trig functions of interest.

[Tyler Josephson ⚛️]

Thanks to Eric Wieser for suggesting that we bundle all our instances together into RealLike.