How to resolve the issue with signed zero?

I came to know about 'signed zeroes' only now as I am trying to deal with complex numbers. Here is the problem:

PROGRAM sZ
       IMPLICIT NONE
       REAL(KIND(0.d0)),PARAMETER :: r=0.2
       COMPLEX(KIND(0.d0)) :: c0
       c0=(0.0,0.5); print*,sqrt(c0**2-r)
       c0=(-0.0,0.5); print*,sqrt(c0**2-r)
       END PROGRAM sZ

The sign of the imaginary part changes.

(0.0000000000000000,0.67082039547127081)
(0.0000000000000000,-0.67082039547127081)

Any instructions/suggestions to resolve this issue are very much welcome.

Solution 1:

The values you see follow from the intermediate step and the specification of the sqrt intrinsic.

The final result depends on the value of the argument to sqrt: as the result has real part the sign of the result's imaginary part is that of the sign of the argument's imaginary part.

What are the signs of the imaginary parts of the two arguments? The imaginary parts of (0.,0.5) and (-0.,0.5) squared are positive in the first case and negative in the second (with the expected implementation of signed zero). Subtracting a real entity does not affect the sign of the imaginary part of the object which becomes the argument to sqrt.

At one level, you are getting the correct result so the only thing to resolve is your expectation of getting the same result. However, consider

  implicit none

  complex :: z1=(-1.,0.), z2=(-1.,-0.)

  print*, z1==z2, SQRT(z1)==SQRT(z2)
end program

It is easy to see why this is confusing even if correct. If we don't want this to happen, we can force any zero imaginary part to be of one particular sign (say, non-negative):

! Force imaginary part away from -0.
if (z%Im==0.) z%Im = 0.

(A comment for this is very necessary.)

We can access and set the imaginary part of a complex variable using the z%Im designator, but this is not possible for general expressions. You may want to create a function for more general use:

  implicit none

  complex :: z1=(-1.,0.), z2=(-1.,-0.)

  print*, z1==z2, SQRT(z1)==SQRT(z2)
  print*, z1==z2, SQRT(design_0imag(z1))==SQRT(design_0imag(z2))

contains

  pure complex function design_0imag(z) result(zp)
    complex, intent(in) :: z
    zp = z
    if (zp%Im==0.) zp%Im = 0.
  end function design_0imag

end program

As shown by other answers, there are several ways to implement logic replacing -0 with 0 (leaving other values unchanged). Remember that you experience this only when the result is purely imaginary.

How to resolve the issue with signed zero?

Solution 1:

Related

Recent Posts