# How can I define unique epsilons with bounded range for floating point analysis?

I'd like to use Mathematica for analyzing error bounds in floating point code. To do this, we can use the standard model for floating point error:

$$flt(a\oplus b) = (a\oplus b)(1+\epsilon_i)$$

Where $$|\epsilon_i| < k$$ where k is some small constant (depending on the operation).

So, I could write FltAdd, FltSub, FltMul, FltDiv, etc functions that compute the result with the extra error term. but the trick is that the epsilon for each separate term has to be unique, so e.g. I can't use $$\epsilon_0$$ for A+B and C+D because in actuality they'll have different errors. If I used the same error, an expression such as (A+B)/(C+D) would have the error cancel which isn't what happens.

So, I think I need to somehow build error terms with unique symbols, but I need them to be predictable enough that I can bound all the "sum" errors by e.g. $$u$$ and all the div terms by e.g. $$u-2u^2$$.

But, I'm not nearly good enough at Mathematica to do that so I thought I'd see what people suggest in terms of approaching this.

• Do you know about interval arithmetic? Oct 30, 2022 at 19:31
• I'm aware of the concept at least, would you go that route and just compute on an interval type instead?
– gct
Oct 30, 2022 at 19:50
• Yes, the rules are already built-in. See Interval Arithmetic Oct 30, 2022 at 19:51
• I'm not sure interval arithmetic is what I want. It seems more suited to tracking the error associated with one evaluation of the function, whereas I want to find a bound for the whole expression for all possible inputs.
– gct
Oct 30, 2022 at 21:06

$Version (* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *) Clear["Global*"]  Use Around (EXPERIMENTAL function introduced in v12.0) Around[a, ϵ]  Format[ϵ[n_]] := Subscript[ϵ, n] {#, #[Around[a, ϵ[1]], Around[b, ϵ[2]]]} & /@ {Plus, Subtract, Times, Divide, Power} // Grid[#, Frame -> All] &  Assuming[ϵ > 0, {#, #[Around[a, ϵ]]} & /@ {Sin, Tan, Sqrt, Exp, Log} // Simplify] // Grid[#, Frame -> All] &  EDIT : From the documentation, "For linear computations, Around[x, δ] behaves like a number whose values are distributed according to the normal distribution NormalDistribution[x, δ]." TransformedDistribution[x + y, {x \[Distributed] NormalDistribution[a, ϵ[1]], y \[Distributed] NormalDistribution[b, ϵ[2]]}] ` • Excellent answer. I think directly used this will overestimate the error, since e.g. for adding two values with (1+$\epsilon$) error on them, I'd expect the error to be bounded by$3\epsilon/(1-3\epsilon)\$ but it gives me a good direction to start looking in.