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

2 次查看(过去 30 天)
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 个评论
Stephen23
Stephen23 2020-6-15
编辑:Stephen23 2020-6-15
"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.

请先登录,再进行评论。

回答(1 个)

Matt J
Matt J 2020-6-14
编辑:Matt J 2020-6-14
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 个评论
James Tursa
James Tursa 2020-6-16
编辑:James Tursa 2020-6-16
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 2020-6-16
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Logical 的更多信息

产品


版本

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by