Verified software floating point

Stream: general

[Geoffrey Irving]

I've formalized software floating point interval arithmetic. Currently, I have ordered field operations (+, -, *, /, <, abs), power series evaluation of exp, log, and powers, and a few extras (scale by a power of two). The main types are

1. Floating: 64+64-bit normalized floating point, wheren : Int64 and s : UInt64 represent the real number n * 2^(s - 2^63). n is normalized (has high bit set) unless s = 0. There are nan values to handle overflow.
2. Interval: Two floating point numbers [lo, hi] representing a nonempty interval (possibly [nan, nan] representing univ).

Operations on Interval are all conservative (guaranteed to approximate exact real computations). Everything is built on top of fixed precision integer operations with no new axioms (there is no trust of the native Float type currently). There are (almost) no results saying that the intervals are tight, so it's empirical whether they expand too much, and in various places I'm not happy with the precision (in particular for power series, though I know some of the reasons why).

I'm going to use this next for two things:

1. An approx tactic which proves real constant formulas by moving them over to Interval.
2. Some Mandelbrot set renders.

I'd be curious if anyone has recommendations for the best way to implement the approx tactic. I could do it by pure pattern matching, but my guess is pushing a bunch of things into a DiscrTree is best. If anyone has a paragraph-length suggestion of how they'd go about a prototype for that, I'd be eager to see get advice.

[Patrick Massot]

Did you look at the existing tactics in Coq and Isabelle?

[Geoffrey Irving]

Currently I have a bunch of lemmas registered for mono which lets me dispatch interval computations. For example,

@[mono] lemma mem_approx_mul [Mul R] [Mul A] [Approx A R] [ApproxMul A R] {a b : R} {x y : A}
    (ax : a ∈ approx x) (yb : b ∈ approx y) : a * b ∈ approx (x * y) := by
  apply subset_approx_mul (singleton_subset_iff.mpr ax) (singleton_subset_iff.mpr yb)
  simp only [mul_singleton, image_singleton, mem_singleton_iff]

[Patrick Massot]

Aren't those eligible to become gcongr lemmas?

[Geoffrey Irving]

I should clarify that I know conceptually how to implement the tactic, so I'm asking if people have advice on Lean details. But admittedly I haven't formulated a very clear question.

[Geoffrey Irving]

No, gcongr is very restrictive: the types have to match on the two sides. Here the congruence is between reals and Interval.

[Geoffrey Irving]

E.g., if I want to turn the addition of two rational numbers into intervals, I need to call Interval.ofRat twice, use its associated approximation lemma Interval.approx_ofRat to show the intervals contain the results, then use approx_add or mem_approx_add to show that the addition is conservative.

[Heather Macbeth]

Have you compared the performance as an aesop wrapper rather than as an apply_rules wrapper (which is what mono is)?

[Heather Macbeth]

Sample code setting it this kind of thing as an aesop wrapper: https://github.com/teorth/pfr/blob/master/PFR/Tactic/Finiteness.lean

[Geoffrey Irving]

Thank you! No, I haven't profiled mono here, and I do expect the aesop version would be faster.

[Sébastien Gouëzel]

Two Isabelle references about their corresponding tactic: https://www.researchgate.net/publication/238740304_Proving_Inequalities_over_Reals_with_Computation_in_IsabelleHOL for the short version, and https://home.in.tum.de/~hoelzl/documents/hoelzl09diplomathesis.pdf for the expanded version.

[Wenda Li]

For Coq interval arithmetic, it may be worth having a look at https://coqinterval.gitlabpages.inria.fr, which, in my impression, has better performance than the Isabelle implementation suggested by Sébastien. Also, Isabelle has affine arithmetic that could be useful when facing linear dependence between variables.

[Jannis Limperg]

The Coq interval tactic clocks in at ~40k LOC (a third of Lean 4), so don't hold your horses.

(Edit: I just realised that this is entirely the wrong metaphor. :upside_down: Well, you know what I mean.)

[Geoffrey Irving]

I mean, I've already done all the hard work: I have fully formalized interval arithmetic. The tactic is just calling the appropriate lemmas.

[Geoffrey Irving]

But probably I am going to get distracted and make a Mandelbrot picture first... :)

[Sébastien Gouëzel]

The reference I gave is among other things interesting for practical tricks to compute quickly sqrt, arctan, cos and so on in a pretty efficient way, better than naive truncation of the Talor series (see Section 4.4 there)

[Geoffrey Irving]

Yes, it's possible to have way tighter intervals than what I'm doing, but if one just wants 10 significant digits the basic stuff is okay. And speed is not important for evaluating simple formulas (it'll be more important in the heavy rendering applications).

[Jannis Limperg]

The Coq tactic (paper) uses reflection with a specialised integer representation and does something called "bisection". But I can't tell how important these features are.

[Geoffrey Irving]

bisection presumably just means chopping the intervals in half and reevaluating. One typically doesn't need to do that if one starts with point-like intervals: it's for proving that large regions are correct. I'm mostly intending to use this in the point-like case, including for rendering, and using theory to give me the region stuff (in particular the Koebe quarter theorem for Mandelbrot renders).

[Geoffrey Irving]

But the code does work sensibly for large intervals. For division, exp, and log I handle them by using monotonicity to reduce to the point-like case, though.

[Geoffrey Irving]

I think reflection just means the normal "proof by reflection" approach, which is also what I'm referring to here.

[Geoffrey Irving]

E.g., to evaluate exp [a,b] I do separate argument reduction and Taylor series evalutions for exp a and exp b (since the multiple of log 2 one wants to subtract may be different) and then take the union.

[Jannis Limperg]

Geoffrey Irving said:

I think reflection just means the normal "proof by reflection" approach, which is also what I'm referring to here.

Yes, it's proof by reflection. This doesn't tend to work so well in Lean because the Lean kernel is very slow compared to Coq. However, if the tactic mostly does computations on Nat, which is builtin and backed by GMP-style bignums, then proof by reflection is maybe fine.

[Geoffrey Irving]

Yes, for the Mandelbrot stuff I'll get around that via by trusting the compiler.