complex numbers multiplication in double precision

Hello,

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

Carlo

Answer 1

[Walter Roberson]

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.

[FastCar]

Ok, I have missed an important detail. The information I have is the phase. So the first equation uses cos and sin to get the real part and the imaginary part. I can tell you that doing a 10^7 addition on the multiplications the two methods give a different result. My idea was that the second one is more precise because does less computations and thus less approximations, but this is far from being a demonstration that it is more accurate. Any idea how to demonstrate that it should be more accurate?

[Walter Roberson]

I need the formulae you used.

Answer 2

[James Tursa]

"... 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 =
    1.000000000000000e+010
>> 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 =
10000000000.0
>> vp = vpa(pi)
vp =
3.1415926535897932384626433832795
>> n = floor(va/(2*vp))
n =
1591549430
>> mod(a,2*pi)
ans =
   5.773954618557402
>> double(va-n*2*vp)
ans =
   5.773954235013852
>> 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.

Answer 3

[FastCar]

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.

[James Tursa]

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.