OT: floats
Quirin Schroll
qs.il.paperinik at gmail.com
Fri Jun 14 12:15:59 UTC 2024
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.
More information about the dip.development
mailing list