CODEWITHSUNDEEP
matlabfloating-pointprecisionfloating-accuracy
Elizabeth Welz
Elizabeth Welz

Reputation: 41

Accuracy of the computation

What will be the computation accuracy if both x and y can be stored (exactly) in floating-point in the system.

x - (x/y) * y

Upvotes: 3

Views: 104

Answers (1)

aka.nice
aka.nice

Reputation: 9362

In case of underflow, the error can be as large as x for example if the ratio x/y completely vanishes and is rounded to zero.

For example, take x=2^-100, y=2^1000.

Otherwise, if x/y is not denormalized (does not underflow), my guess is that you will get a perfect zero most of the time, and at most eps(x) - that is 1 ulp(x) from time to time, as long as Matlab engine respects IEEE754 standards.

The reason is this one: let's note the rounding error e

float(x/y) = x/y + e

With round to nearest, tie to even IEEE754 mode, we have

abs(e) < ulp(x/y)/2

For division, this is a strict inequality, otherwise the division would have to be exact. And if it is exact, we are sure that the quotient fits in available significand (unless the division underflows).

When we multiply back by y, we get this exact result:

x + e*y

We need to have e*y >= ulp(x)/2 for this exact result to be rounded off x
That is y*ulp(x/y)/2 > e*y >= ulp(x)/2.
This can happen, for example try:

x=2^52+2^51+1 , y=2^52+1 , x - (x/y)*y , y*eps(x/y)/2 > eps(x)/2

But y*ulp(x/y)/2 > e*y >= 3*ulp(x)/2 cannot, so the result cannot be 2 ulp off.
This is tedious to demonstrate, so my answer will remain qualified as a guess.

The last subtraction operation will be exact - see https://en.wikipedia.org/wiki/Sterbenz_lemma

Upvotes: 5

Related Questions