MATLAB Answers

Why is a Number Divided by Itself Not Equal to Unity?

3 views (last 30 days)
Paul
Paul on 14 Jun 2020
Commented: Paul on 16 Jun 2020
w=0:(pi/100):pi;
w = w([17 66 93]);
x = exp(1i*w);
x./x
ans =
1.0000e+00 + 4.8645e-17i 1.0000e+00 - 4.9461e-17i 1.0000e+00 - 2.6884e-17i
I imagine this result has something to do with how Matlab implements division of complex numbers and I realize that the error is less than eps(1), But I thought it was an interesting result and it would never have occurred to me that this could happen. I checked the doc page for rdivide, but it doesn't say anything about its algorithm when both arguments are complex. Just curious if anyone can explain why this result happens (it only happened for these three particular values of the original w vector).

  6 Comments

Show 3 older comments
Paul
Paul on 14 Jun 2020
Implementing the complex division exactly as on that matworld page gets the expected answer:
>> a=real(x);b=imag(x);c=real(x);d=imag(x);
>> (a.*c+b.*d)./(c.^2+d.^2) +1j*(b.*c - a.*d)./(c.^2+d.^2) - 1
ans =
0 0 0
So I'm not sure what side effect you're referring to.
Stephen Cobeldick
Stephen Cobeldick on 15 Jun 2020
"Implementing the complex division exactly as on that matworld page gets the expected answer:"
Sure, but I was not making any comment on that specific algorithm: it is quite possible that the actual implentation (in MATLAB, BLAS, or whatever) would be something algebraically equivalent but that uses fewer/faster/more stable/whatever operations, and hence has different accumulated floating point error to what you tried.
My point was simply that complex division is not a single atomic numeric operation, so whatever algorithm that is used to implement it could quite conceievably accumulate some floating point error: the authors of numeric libraries must always compromise speed vs memory vs accumulated error vs stability vs whatever else.

Sign in to comment.

Answers (1)

Matt J
Matt J on 14 Jun 2020
Edited: Matt J on 14 Jun 2020
It depends what algorithm is used to perform complex division x./y. For example, if it is implemented as x./abs(y).^2.*conj(y), then there will be floating point errors introduced in these operations, even when x=y,
>> x./abs(x).^2.*conj(x) - 1
ans =
1.0e-15 *
0 0 0.2220
If it had been done, however, by conversion to complex exponentials, I don't think we should have seen residual errors,
>> [r,t]=deal(abs(x),angle(x));
>> (r./r).*exp(1i*(t-t))-1
ans =
0 0 0

  10 Comments

Show 7 older comments
Paul
Paul on 16 Jun 2020
I don’t know much, or anything really, about language standards, but in the little bit of reading I did about the FORTRAN standard I didn’t see any requirement that complex division must yield 1 for x/x. Does any other language have that requirement in a standard? At the risk of asking a stupid question, do you know if any language has a requirement that x/x must equal 1 for non-complex types? Or is that a requirement of the IEEE 754 standard?
James Tursa
James Tursa on 16 Jun 2020
These are not supid questions IMO.
You would need to look at IEEE floating point standards for stuff like x/x = 1. I think IEEE requires that the calculation result be correct as if the calculation is done in infinite precision and rounded according to the lsb rounding rule in effect (e.g., typically "round to even"), so that would mean x/x = 1 would be required to be true for all possible real finite floating point bit patterns (except 0). It would not be hard to write a MATLAB program to test every possible single and double bit pattern to see if it is true. (But I would be shocked to find out that any such bit pattern failed this test ... this is implemented at the chip level)
For language standards regarding complex arithmetic, I am unaware of any particular requirements (e.g., z/z = 1 or avoiding overflow & underflow when possible). I think that Java has some floating point replicability requirements in their standards but I don't know the details.
Paul
Paul on 16 Jun 2020
Thank you both for the discussion. I did some research and learned about three algorithms for complex division, among other things. I'll probably never have to code them myself, but still good to know they exist.

Sign in to comment.

Products


Release

R2019a

Translated by