sqrt problem... or compiler?

[Oleg]

I have problem with calculation of sqrt of complex argument. Can anybody explain whats wrong.

there is my simple code for calculation of sqrt( (-100,0) )

program tupost
complex, parameter :: J = (0,1)
complex a, b, ax

ax = (10.0, 0.0) ! I need this paarmeter to be complex
a = (-ax**2)

b = sqrt(real(a) + J*imag(a))
print *,' b = ', b

!b = sqrt(a)
!print *,' b = ', b

end program tupost

This gives correct result b = (0,10)

However, direct calculation b = sqrt(a) gives (0,-10). Moreover, if you remove comments for calculation together sqrt(a)

you get (0,-10) for this and previous calculation!!! Is it correct, or its a compiler bug?

I use VFC version 10.1.025

[ender01]

Quoting - magamp

I have problem with calculation of sqrt of complex argument. Can anybody explain whats wrong.

there is my simple code for calculation of sqrt( (-100,0) )

program tupost
complex, parameter :: J = (0,1)
complex a, b, ax

ax = (10.0, 0.0) ! I need this paarmeter to be complex
a = (-ax**2)

b = sqrt(real(a) + J*imag(a))
print *,' b = ', b

!b = sqrt(a)
!print *,' b = ', b

end program tupost

This gives correct result b = (0,10)

However, direct calculation b = sqrt(a) gives (0,-10). Moreover, if you remove comments for calculation together sqrt(a)

you get (0,-10) for this and previous calculation!!! Is it correct, or its a compiler bug?

I use VFC version 10.1.025

Hi Magamp,

Part of the 2003 standard includes changes to intrinsic functions atan2, log, and sqrt for processors that distinguish between positive and negative real zeros (pretty much everything). In particular, for complex arguments now returns a negative imaginary result if the real part of the result is zero and the imaginary part of the argument is less than zero. That being said, at first glance, the imaginary part of your argument is not less than zero. So I looked at the binary representation of a:

11000010110010000000000000000000
10000000000000000000000000000000

Under the IEEE 32 bit standard, the leading bit is the sign bit and both numbers are negative. The sign flip on the imaginary part is probably due to the sign change forced by the unary operator "-":

ax**2:

1000010110010000000000000000000
0

-ax**2

11000010110010000000000000000000
10000000000000000000000000000000

When you construct the argument using the parameter, the argument becomes:

11000010110010000000000000000000
0

So no, it's not a compiler bug.

Cheers,

Rich

[Oleg]

Thank you,

Is it normal that

b = sqrt(real(a) + J*imag(a))
print *,' b = ', b

gives (0,10)

but two lines together

b = sqrt(real(a) + J*imag(a))
print *,' b = ', b
b = sqrt(a)
print *,' b = ', b

result

(0,-10)

(0,-10)

??

why in this case the imaginary part of argument using the parameter is changed?

[ender01]

Quoting - magamp

Thank you,

Is it normal that

b = sqrt(real(a) + J*imag(a))
print *,' b = ', b

gives (0,10)

but two lines together

b = sqrt(real(a) + J*imag(a))
print *,' b = ', b
b = sqrt(a)
print *,' b = ', b

result

(0,-10)

(0,-10)

??

why in this case the imaginary part of argument using the parameter is changed?

Well, yes and no. There's probably an optimizer that is essentially rewriting your code so what you've programmed in may not exactly be the exact thing that's being executed. For instance, the optimizer will rearrange the code to do things like increase execution speed or increase precision. In the latter case, the optimizer will try to keep the results of sequential computations in a double length register for instance so that you'll see less round-off error. Print statements in the middle of a sequence of computations, however, will force the value to be stored in it's declared precision so you'll see differences with and without print statements in some cases due to increased round-off. That said, I can't say for certain that that's what's going on but it does fit (though in this case I'd suspect code selection, variable simplification rather than inlining). When I compile in with the debug settings (optimization off) I don't see that behavior but when I compile with the release settings (optimization at a higher level ) I do see it so I suspect the optimizer. If it's a problem, you can lower the optimization setting until you get the right result.

Cheers,

Rich

[Oleg]

Thank you. Now it is clear.