Casting to and from Floats

Stream: new members

[Caleb Schultz Kisby]

Hi again,

I'd like to give my program Floats, and then maybe convert them to something like Rat or Real that's easier to reason about. I don't mind if the things I prove about Rat or Real aren't actually true of the Float, it's an illusion I'm happy to maintain. (At the end of the day, I'd like to be able to both evaluate my Floats and prove things about numbers that they approximate.)

I've looked around in the Zulip archive, and I've found several threads where others mention using Reals instead of Floats. But I'm having trouble figuring out how to actually cast between them. I've looked in Init.Data.OfScientific, Mathlib.Data.Real.Basic, and Std.Data.Rat.Basic.

Am I just looking in the wrong places? Am I approaching this the wrong way?
Thanks,
Caleb

[Tomas Skrivan]

To achieve this I currently take an axiomatic approach. I just create a new type that wraps around Float and state that it is a field. Here is a code that just postulates that such type is commutative ring

import Mathlib.Algebra.Ring.Basic

structure MyReal : Type where
  val : Float

def _root_.Int.toFloat : Int → Float
| .ofNat z' => z'.toFloat
| .negSucc z' => -(z'+1).toFloat

instance commRing : CommRing MyReal := by
  refine' { natCast := fun n => ⟨n.toFloat⟩
            intCast := fun z => ⟨z.toFloat⟩
            zero := ⟨(0 : Float)⟩
            one := ⟨(1 : Float)⟩
            mul := λ x y => ⟨x.1 * y.1⟩
            add := λ x y => ⟨x.1 + y.1⟩
            neg := λ x => ⟨-x.1⟩
            sub := λ x y => ⟨x.1 - y.1⟩
            npow := λ n x => ⟨x.1^n.toFloat⟩
            nsmul := λ n x => ⟨n.toFloat * x.1⟩
            zsmul := λ z x => ⟨z.toFloat * x.1⟩,
            .. }
  all_goals admit  -- Just assume that all the laws are true

instance : OfScientific MyReal := ⟨λ m s e => ⟨(OfScientific.ofScientific m s e)⟩⟩
instance : ToString MyReal := ⟨λ x => x.1.toString⟩

#eval (42.0 : MyReal) + (0.123456 : MyReal)

instance : LE MyReal := ⟨λ x y => x.1 ≤ y.1⟩
instance : LT MyReal := ⟨λ x y => x.1 < y.1⟩
instance : DecidableRel (λ x y : MyReal => x ≤ y) :=
  λ x y => if h : x.1 ≤ y.1 then isTrue h else isFalse h
instance : DecidableRel (λ x y : MyReal => x < y) :=
  λ x y => if h : x.1 < y.1 then isTrue h else isFalse h

example
  : (1.0 : MyReal) + (-1.0 + 0.00000000000000001 : MyReal) ≤ 0.0 := by native_decide

example
  : (1.0 : MyReal) + (-1.0 + 0.00000000000000001 : MyReal) > 0.0 := by rw[← add_assoc]; native_decide

The end shows that you can get into trouble once you start playing around with native_decide and decidable inequalities. Here is a thread discussing this inconsistency.

[Tomas Skrivan]

Another approach I'm experimenting with now, is where I provide the compiler with different definitions of operations on reals. This can be done with impelmented_by attribute.

import Mathlib.Data.Real.Basic
import Mathlib.Data.Complex.Exponential

def MyReal := Real

-- `MyReal` and mathlib's reals `ℝ` are defeq
example : MyReal = ℝ := by rfl

-- To convert from float to real we just perform unsafe cast
-- We have to make sure that any operation on `MyReal` has its implementaion
-- overridden by corresponding implementation with Floats
unsafe def MyReal.toFloat (x : MyReal) : Float := unsafeCast x
unsafe def Float.toMyReal (x : Float)  : MyReal := unsafeCast x

namespace MyReal

unsafe def Impl.add (x y : MyReal) : MyReal := (x.toFloat + y.toFloat).toMyReal

@[implemented_by Impl.add]
def add (x y : MyReal) : MyReal := ⟨x.1 + y.1⟩

unsafe def Impl.ofScientific (mantissa : Nat) (exponentSign : Bool) (decimalExponent : Nat) : MyReal :=
  (OfScientific.ofScientific mantissa exponentSign decimalExponent : Float).toMyReal

instance : Add MyReal := ⟨add⟩

@[implemented_by Impl.ofScientific]
def ofScientific (mantissa : Nat) (exponentSign : Bool) (decimalExponent : Nat) : MyReal :=
  let exp := (10:ℝ)^decimalExponent
  if exponentSign then
    (mantissa / exp : ℝ)
  else
    (mantissa * exp : ℝ)

instance : OfScientific MyReal := ⟨ofScientific
⟩
unsafe instance : ToString MyReal := ⟨λ x => x.toFloat.toString⟩

#eval (42.0 : MyReal) + (0.123456 : MyReal)

However, I'm unable to show that MyReal is a field without doing lots of work.

[Caleb Schultz Kisby]

Ah, thank you so much! This is perfect. The inconsistency with native_decide is fine, and is what I'd expect (since you're evaluating the Floats). I'll just avoid evaluating Float in the proofs.

[Caleb Schultz Kisby]

Also, I didn't know about native_decide until now (I'm back to learning Lean again!). I used to use Agda, and so I'm unfamiliar with the tactic language. But things like this make it hard to go back to just terms :-)

[Kevin Buzzard]

If you don't want to evaluate the floats then why don't you just use reals instead?

[Tomas Skrivan]

Evaluating floats via native_decide in a proof and running a program with floats are two different things.

[Caleb Schultz Kisby]

If I can ask a follow-up question -- is there any official plan for handling this sort of cast/conversion in the Lean4 mathlib? (Again, I'm new to Lean, so I'm not familiar with mathlib development. I'm just curious, and I'd also like to avoid writing code that will look archaic in 5 years!)

[Martin Dvořák]

Mathematicians are not interested in floats afaik. Some library other than Mathlib4 might incorporate it tho.

[Caleb Schultz Kisby]

Oh whoops, good point. I forgot that Float is dealt with in Init.Data.Float and Init.Data.OfScientific anyway :sweat_smile:

[Kevin Buzzard]

For sure mathlib is not interested in defs which make the system inconsistent :-)

[Reid Barton]

Applied mathematicians of a certain type are plenty interested in floats. If we come up with some useful way to talk about the relationship between idealized calculations with real numbers and actual algorithms with floats, I think it would be within scope for mathlib. But of course it would have to be done in a way that doesn't involve adding false axioms.
At the moment there isn't any plan to do so though, or even any concrete idea of how to do so really.

[Kevin Buzzard]

Presumably the way to do this is to have some kind of interval arithmetic and a map from Float to that? Did anyone ever write down precise functions lower_bound and upper_bound from Float to Rat which have some kind of useful meaning?

[Eric Wieser]

Do we have docs4#Float.toRat ?

[Eric Wieser]

I'm not sure the lower/upper bound is that interesting, it's probably just half-way between the adjacent Floats, perhaps with a different definition for each rounding mode

[Eric Wieser]

Is it even possible to prove things about Float in lean4? It looks like you can't actually access the underlying bits

[Mario Carneiro]

you can prove things about concrete Floats using native_decide (which can see through opaque and hence is a proper axiomatic extension), but that's not good enough for doing general properties about float representation. I would probably do full IEEE modeling by having a FP.Float type a la docs4#FP.Float and axiomatize the bijection between them (this may require additional functions in core if necessary to define the bijection)

[Eric Wieser]

Why would we need the bijection to be axiomatized? We don't have any such axiomatization for things like Uint64

[Mario Carneiro]

because Float is opaque

[Mario Carneiro]

I'm talking under the assumption that the definition of Float in core does not change

[Mario Carneiro]

The type docs4#FloatSpec which gives the properties of the float type itself are not sufficient to establish e.g. its cardinality

[Mario Carneiro]

you can introduce the "to" and "from" functions using opaque, similar to most float ops, but any assertions about those functions would need to be axiom, for example (Float.add x y).toFP = FP.Float.add x.toFP y.toFP

[Mario Carneiro]

It would be possible to have structure Float where ofFP :: toFP : FP.Float though in core; the toFP and ofFP functions would be overridden in this case such that Float has the ABI of a f32 while FP.Float is a regular lean inductive containing the "unpacked" components

[Eric Wieser]

I think that was the version I was expecting

[Eric Wieser]

But maybe that's asking too much of core especially if we want implementations of all the float functions within the theory

[Martin Dvořák]

This is a long shot, but...

I'd like to do computation with matrices over approximate complex numbers (where the real part is a Float and the imaginary part is a Float as well) and I would prefer not to define everything myself. Does Lean 3 or Lean 4 provide anything useful for my purpose?

[Ruben Van de Velde]

I think the short answer is "no" :)

[Ruben Van de Velde]

Not sure if there's a longer answer yet

[Eric Wieser]

I've argued before that we should have such a construction, since it would absorb docs#gaussian_int

[Eric Wieser]

In the meantime, you could use a clifford_algebra!

[Reid Barton]

You should check out https://github.com/lecopivo/SciLean

[Henrik Böving]

Martin Dvořák said:

This is a long shot, but...

I'd like to do computation with matrices over approximate complex numbers (where the real part is a Float and the imaginary part is a Float as well) and I would prefer not to define everything myself. Does Lean 3 or Lean 4 provide anything useful for my purpose?

If you dont want to prove anything you can just do Array (Array (Prod Float Float)) orr am I overlooking something?

[Reid Barton]

I think this is "defining everything myself"

[Eric Wieser]

I would imagine that you can use docs4#Matrix even if you can't use docs4#complex. Though for computation you'll want to occasionally un-lazy the matrix by converting to nested arrays and back

[Reid Barton]

Probably closer to "always" than "occasionally".

[Reid Barton]

The Mul instance needs some well-behavedness of the underlying ring, which you would need to assert (even though it's false).

[Reid Barton]

If you don't want to prove theorems then rolling your own seems better in pretty much every respect

[Eric Wieser]

Reid Barton said:

The Mul instance needs some well-behavedness of the underlying ring, which you would need to assert (even though it's false).

Ah, right; because we use finset.univ.sum to collect the terms, but that's only well-defined if addition is associative.

[Reid Barton]

Which is unavoidable I think, because it only assumes a fintype instance on the indexing types

[Reid Barton]

Also, some good old-fashioned for loops are likely to be a lot more efficient than finset.univ.sum, depending on how good Lean is at inlining

[Reid Barton]

When the body of the loop is two fetches and an fmadd, there is a lot of room to make it slow.