complex numbers multiplication in double precision

108 views (last 30 days)
FastCar on 1 Jul 2017
Commented: James Tursa on 3 Jul 2017
I have to multiply couple of complex numbers and then I have to add all the product.
Let's say the two complex numbers (with unitary module) are c1 and c2. What I need is prod=c1 x c2* (but I can substitute c2 with its conjiugate in our example). c1 has phase phi1, c2 has phase phi2. I can do it in chartesian coordinates
c1 x c2* = real(c1) x real (c2) + imag(c1) x imag(c2) + j x (imag(c1) x real (c2) - real(c1) x imag(c2)
(if I have used the right formula)
or I can use
c1 x c2* = cos(phi1-phi2) + j x sin(phi1-phi2)
Mind that the phases are of the order of magnitude of 10^10 radians.
The second way is of course faster, but my question is which one is the more precise given that the computation is done in double precision?
Mind also that I have to sum 10^7 products and I can tell you that there is a difference, a not neglectable difference.
Many thanks in advance

Answers (3)

Walter Roberson
Walter Roberson on 1 Jul 2017
Edited: Walter Roberson on 1 Jul 2017
On my system, even excluding the time it takes to calculate the angle, the cos/sin version is about 1.8 to 2.4 times slower than the other version. This is to be expected as cos and sin require a many more calculations. There are hardware instructions, yes, but they are slower than plain multiply or addition.
The second version appears to have marginally better numerical properties eps(1) compared to 2*eps(1). That is not really a meaningful difference compared to the rough-off error in the rest of the calculations.

Sign in to comment.

James Tursa
James Tursa on 2 Jul 2017
Edited: James Tursa on 2 Jul 2017
"... Mind that the phases are of the order of magnitude of 10^10 radians. ..."
When you have angles that large, I have to start questioning where they came from and how exactly you are using them in your code. Realize that 10^10 is more than 9 orders of magnitude larger than 2pi, so this will come into play in the "accuracy" you expect for calculations. E.g., a naive example:
>> a = 10^10
a =
>> cos(a)+sin(a)*i
ans =
0.873119622676856 - 0.487506025087511i
>> cos(mod(a,2*pi))+sin(mod(a,2*pi))*i
ans =
0.873119809656583 - 0.487505690208075i
So you can see a difference in the 6th digit. What is going on? It turns out that the trig functions have a very sophisticated method of dealing with the argument. I don't know the exact method that is used behind the scenes, but effectively it treats the IEEE double precision argument as if it has the same absolute precision as a IEEE double precision pi. E.g.,
>> va = vpa(a)
va =
>> vp = vpa(pi)
vp =
>> n = floor(va/(2*vp))
n =
>> mod(a,2*pi)
ans =
>> double(va-n*2*vp)
ans =
>> cos(ans)+sin(ans)*i
ans =
0.873119622676856 - 0.487506025087511i
So you can see by treating the angle as if it had the same absolute precision as a double precision pi value we can do the mod in such a way that we can uncover how cos() and sin() are behaving behind the scenes. It is effectively a different angle that is used compared to the naive mod(a,2*pi) angle result.
Bottom line is that with 10^10 radian angles you need to be very careful that your problem has real meaning at these values, and then be extra careful how you are doing calculations with them to come up with meaningful results. For your phi1-phi2 method, realize that the cos() and sin() functions are going to treat that phi1-phi2 double precision result as if it had absolute precision the same as a double precision pi. This effect may be creeping into your results.

FastCar on 2 Jul 2017
Edited: FastCar on 2 Jul 2017
I try to clarify my question: the phase is given by the distance of a satellite from the Earth multiplied by
4 x pi / lambda
The distance is big (around 40000 km) and thus gives a huge phase. So far I used these expressions:
distance = norm ( position vector )
phi = ( 4 x pi x distance ) / lambda
method one
c1= cos(phi1) + j x sin(phi1)
c2= cos(phi2) + j x sin(phi2)
c1 x c2* = cos(phi1) x cos(phi2) + sin(phi1) x sin(phi2) + j x (sin(phi1) x cos (phi2) - cos(phi2) x sin(c2)
method 2
c1 x c2* = cos(phi1-phi2) + j x sin(phi1-phi2)
For what I know every computation in matlab should be done in double precision and without the symbolic packed pi should have only 15 figures (don't blame me if I am wrong). Mind that 15 figures (and thus 14 sure figures) is a good output in my point of view.
I hope I have clarified the question.
I haven't answered to the single question because I think both of you needed the full picture.
  1 Comment
James Tursa
James Tursa on 3 Jul 2017
I don't quite understand your physical application yet, but keep in mind what I have already written in my Answer. Because phi is so large, you have effectively already lost 9 significant digits from your answer. So the c1 x c2* result will only have about 6 meaningful digits of accuracy at most.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!