then we should change bit(s) because it means dick/penis in french /s

A min-max-median heap, I would assume.

So would I!

No, I haven't. I have used them with shallow embedding techniques that are relatively similar, but not in this way.

There is nothing inherently difficult about practical implementations of continuous numbers for automated reasoning compared to more discrete mathematical structures. They are handleable by standard FOL itself.

See ACL2's support for floating point arithmetic.

https://www.cs.utexas.edu/~moore/publications/double-float.p...

SMT solvers also support real number theories:

https://shemesh.larc.nasa.gov/fm/papers/nfm2019-draft.pdf

Z3 also supports real theories:

https://smt-lib.org/theories-Reals.shtml

I'm curious about how you'd do that, too. I haven't tried doing anything like that, but I'd think you'd start by trying to formalize the usual optimality proof for the Kalman filter, transfer it to actual program logic on the assumption that the program manipulates real numbers, try to get the proof to work on floating-point numbers, and finally extract the program from the proof to run it.

https://youtu.be/_LjN3UclYzU has a different attempt to formalize Kalman filters which I think we can all agree was not a successful formalization.

The link was to show Coq's floating point library.

> As a concrete example, how would you prove an implementation of a Kalman filter is correct?

This would be fun enough to implement if I had more time. But I need to spend time with my family. I'm usually loathe to just post LLM output, but in this case, since you're looking just for the gist, and I don't have the time to write it out in detail, here's some LLM output: https://pastebin.com/0ZU51vN0

Based on a quick skim, I can vouch for the overall approach. There's some minor errors/things I would do differently, and I'm extremely doubtful it typechecks (I didn't use a coding agent, just a chat window), but it captures the overall idea and gives some of the concrete syntax. I would also probably have a third implementation of the 1-d Kalman filter and would prefer some more theorems around the correctness of the Kalman filter itself (in addition to numerical stability). And of course a lot of the theorems are just left as proof stubs at the moment (but those could be filled in given enough time).

But it's enough to demonstrate the overall outline of what it would look like.

The overall approach I would have with floating point algorithms is exemplified by the following pseudo-code for a simple single argument function. We will first need some floating point machinery that establishes the relationship between floating point numbers and real numbers.

  # Calculates the corresponding real number from the sign, mantissa, and exponent
  floatToR : Float -> R
  
  # Logical predicate that holds if and only if x can be represented exactly in floating point
  IsFloat(x: R)

Then define a set of three functions, `f0` which is on the idealized reals, `f1` which is on floating point numbers as mechanically defined via binary mantissa, exponent, and sign, and `f2` which is on the finite subset of the reals which have an exact floating point representation.

  f0 : R -> R
  f1 : Float -> Float
  f2 : (x : R) -> IsFloat(x) -> R

The idea here is that `f0` is our "ideal" algorithm. `f1` is the algorithm implemented with exact adherence to low-level floating point operations. And `f2` is the ultimate algorithm we'll extract out to runnable code (because it might be extremely tedious, messy, and horrendous for runtime performance, to directly tie `f1`'s explicit representation of mantissa, exponent, and sign, which is likely some record structure with various pointers, to the machine 32-bit float).

Then I prove the exact correspondence between `f1` and `f2`.

  f1_equal_to_f2 : forall (x: R), (isFloatProof: IsFloat(x)) -> f1(floatToR(x)) = floatToR(f2(x, isFloatProof))

This gives the following corollary FWIW.

  f2_always_returns_floats : forall (x: R), IsFloat(x) -> IsFloat(f2(x))

This lets me throw away `f1` and work directly with `f2`, which lets me more easily relate `f2` and `f0` since they both work on reals and I don't need to constantly convert between `Float` and `R`.

So then I prove numerical stability on `f2` relative to `f0`. I start with a Wilkinson style backwards error analysis.

  f2_is_backwards_stable : forall (x: R), exists (y: R), IsFloat(x) -> f2(x) = f0(y) and IsCloseTo(y, x)

Then, if applicable, I would prove a forward error analysis (which requires bounding the input and then showing that overflow doesn't happen and optionally underflow doesn't happen)

  f2_is_forwards_stable : forall (x: R), IsFloat(x) -> IsReasonablyBounded(x) -> IsCloseTo(f2(x), f0(x))

And then given both of those I might then go on to prove some conditioning properties.

Then finally I extract out `f2` to runnable code (and make sure that compilation process after that has all the right floating point flags, e.g. taking into consideration FMA or just outright disabling FMA if we didn't include it in our verification, making sure various fast math modes are turned off, etc.).

This is potentially horrendous for performance, and even worse, unpredictably so. Instead of getting the incorrect answer (or NaN) with floating point if you have an unfortunate series of calculations, you get extreme memory blowup, where your numerator and denominator can explode in size, which in turn leads to runtime slowdown. In the worst case you can actually run out of memory (because certain fairly natural calculations can cause your numerators and denominators to very rapidly explode in size).

Committing to unbounded rationals basically opens your system up to DoS attacks.

You could decide to bound rationals in the numerator and denominator, but then you've more or less reinvented a form of floating point.

In order of preference for a high reliability production system I would use:

1. Integers

2. Fixed point

3. Floating point

4. Rationals [waaaayyyy down in fourth place]

Most systems don't have high enough reliability needs though that I would favor fixed point over floating point and so in practice I rarely use fixed point (simply because library and language support is far worse).

The standardized operations, e.g. multiplication or square root extraction, are precisely defined, i.e. the result is always defined exactly, by the combination of the corresponding operation with real numbers and by the rounding rule that is applied.

IEEE 754 also contains a list of operations that are recommended, but not defined by the standard, such as the exponential function and other functions where it is difficult to round exactly the result.

For the standardized operations, all the possible errors are precisely defined and they must either generate an appropriate exception or produce as result a special value that encodes the kind of error, depending on how the programmer configures the processor.

The standard is perfectly fine. The support of the standard in the popular programming languages is frequently inconvenient or only partial or even absent. For instance it may be impossible to choose to handle the errors by separate exception handlers and it may be impossible to unmask some of the exceptions that are masked by default. Or you may lack the means to control the rounding mode or to choose when to use FMA operations and when to use separate multiplications.

If you enable all the possible exceptions, including that for inexact results, the value of an expression computed with IEEE 754 operations is the same as if it were computed with real numbers, so you do not need to prove anything extra about it.

However this is seldom helpful, because most operations with FP numbers produce inexact results. If you mask only the exception for inexact results, the active rounding rule will be applied after any operation that produces an inexact result.

Then the expression where you replace the real numbers with FP numbers is equivalent with a more complex expression with real numbers that contains rounding operations besides the explicit operations.

Then you have to prove whatever properties are of interest for you when using the more complex expression, which includes rounding operations.

The main advantage of the IEEE 754 standard in comparison with the pathetic way of implementing FP operations before this standard, is that the rounding operations are defined exactly, so you can use them in a formal proof.

Before this standard, most computer makers rounded the results in whatever way happened to be cheaper to implement and there were no guarantees about which will be the result of an operation after rounding, so it was impossible to prove anything about FP expressions computed in such computers.

If you want to prove something about the computation of an expression when more exceptions are masked, not only the inexact result exception, that becomes more complex. When a CPU allows a non-standard handling of the masked exceptions, like flush-to-zero on underflow, that can break any proof.

Just to clarify for others because it’s a tiny bit clickbaity:

The Patriot missile didn’t kill 28 people accidentally, it simply failed to intercept an enemy missile.

And it wasn’t launched on an incorrect trajectory either, the radar was looking at a slightly wrong distance window and lost track. Furthermore, the error only starts having an effect after 100 hours of operation, and it seems to have only been problematic with the faster missiles in Iraq that the system wasn’t designed for. They rushed the update and they did actually write a function to deal with this exact numerical issue, but during the refactor they missed one place where it should have been used.

28 lives are obviously significant, but just to note that there are many mitigating factors.

I have a postit note idea that says simply "typesafe macro assembler".

I've not fleshed this out yet, but I think a relatively simple system would help deal with all the issues you mention in the first paragraph while allowing escape hatches.