OT: floats

[Quirin Schroll]

On Friday, 14 June 2024 at 06:39:59 UTC, Dom DiSc wrote:
> On Thursday, 13 June 2024 at 19:50:34 UTC, Quirin Schroll wrote:
>> [...] If you implement an algorithm, you have to take into 
>> account how rounding errors propagate through the 
>> calculations. The issue is that you can’t do that intuitively. 
>> You just can’t.
>
> Best is to use an interval-type, that says the result is 
> between the lower and the upper bound. This contains both the 
> error and the accuracy information.
> e.g. sqrt(2) may be ]1.414213, 1.414214[, so you know the exact 
> result is in this interval and you can't trust digits after the 
> 6th. If now you square this, the result is ]1.9999984, 
> 2.0000012[, which for sure contains the correct answer.
> The comparison operators on these types are ⊂, ⊃, ⊆ and ⊇ - 
> which of course can not only be true or false, but also ND (not 
> defined) if the two intervals intersect.
> so checks would ask something like "if(sqrt(2)^^2 ⊂ ]1.999, 
> 2.001[)", so they don't ask only for the value but also for the 
> allowed error.

[Note: The following contain fractions which only render nice on 
some fonts, _excluding_ the Dlang forum font.]

Interval arithmetic (even on precise numbers such as rationals) 
fails when the operations blow up the bound to infinity. This can 
happen because of two reasons: The actual bounds blow up (i.e. 
the algorithm isn’t well-conditioned, evidenced by analysis), or 
analysis shows the bounds are actually bounded, but interval 
arithmetic doesn’t see it. An example for the latter case:

Let *x* ∈ [1⁄4, 3⁄4]. Set
 *y* = 1⁄2 − *x* ⋅ (1 − *x*) and
 *z* = 1 / *y*

Then, standard interval arithmetic gives
 *y* ∈ [1⁄2, 1⁄2] − [1⁄16, 9⁄16] = [−1⁄16, 7⁄16]
 *z* ∈ [−∞, −16] ∪ [16⁄7, +∞]

However, rather simple analysis shows:
 *y* ∈ [1⁄2, 1⁄2] − [3⁄16, 1⁄4] = [1⁄4, 5⁄16]
 *z* ∈ [16⁄5, 4]

Of course, interval arithmetic isn’t mathematically wrong, after 
all
 [16⁄5, 4] ⊂ [−∞, −16] ∪ [16⁄7, +∞],
but the left-hand side interval is a useful bound (it’s optimal, 
in fact), but the right-hand side one is near-useless as it 
contains infinity.

The big-picture reason is that interval arithmetic doesn’t work 
well when a value interacts with itself: It counts all 
occurrences as independent values. In the example above, in *x* ⋅ 
(1 − *x*), the two factors aren’t independent. And as you can 
see, for the upper bound of *y*, that’s not even an issue: If 
interval arithmetic determined that *y* is at most 7⁄16 (instead 
of 5⁄16) and therefore *z* is at least 16⁄7 (instead of 16⁄5), 
that’s not too bad. The issue is the lower bound.

Essentially, what one has to do is make every expression 
containing a variable more than once into its own (primitive) 
operation, i.e. define *f*(*x*) = *x* ⋅ (1 − *x*), then determine 
how *f* acts on intervals, and use that. An interval type 
(implemented in D or any other programming language) probably 
can’t do easily automate that. I could maybe see a way using 
expression templates or something like that where it could 
determine that two variables are in fact the same variable, but 
that is a lot of work.