 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 flips.

(0.0000000000000000,0.67082039547127081)
(0.0000000000000000,-0.67082039547127081)

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

2 Likes

Anton Shterenlikht has a paper on this https://arxiv.org/abs/1712.10230, I believe.

Indeed, it very much depends on the quality of the implementation of the runtime libraries whether the signed zero is respected or not. Anton’s article illustrates that nicely. (Tangential to this: see the post by Beliavsky - Fortran code snippets - #22 by Beliavsky - where he illustrates the use of complex numbers to determine the derivative of a real (!) function with just one evaluation.)

This is a really cool technique, although the more robust version of it is forward mode autodiff (using dual numbers) which can handle a significantly larger set of programs.

@gyrtle: Interesting post!

It depends on your problem. If the real part must be positive due to some physical restrictions, the minus sign is just the result of a rounding error. Then simply enforce a positive number by using abs. If negative values are allowed, you need to cope with the situation anyways: The behavior you have observed exists also for tiny but finite real parts, so there is nothing special about the negative zero.

If you explicitly compare whether the real part is zero, take floating point inaccuracies into account. I found the following link quite useful: Comparing Floating Point Numbers, 2012 Edition | Random ASCII – tech blog of Bruce Dawson. My implementations are on DAMASK/prec.f90 at master · eisenforschung/DAMASK · GitHub. For example

c0 = merge(cmplx(0.0,0.5), cmplx(a,0.5), dEq0(a))

If you test whether the real part is smaller than zero, take into account that -0.0<0.0 gives you .false.. The sign function can help:

real :: n = -0.0, p = +0.0

print*, sign(1.0,n)<0.0, n<0.0
print*, sign(1.0,p)<0.0, p<0.0

If you’re going to post the exact same question on stackoverflow, why not wait for an answers.

I guess what OP want is a one line solution to ensure the result of sqrt(c0**2-r) have positive imaginary part.
Perhaps someone can have such a one line solution?

Of not, OP can simply define a function called positive_imagarypart_sqrt or something, to ensure that

positive_imagarypart_sqrt( xxxxxx )

gives a positive imaginary part.

By the way welcome!

1 Like

Thanks. The paper is very helpful.

1 Like

Sorry for cross-posting. I found this group after I posted on Stackoverflow.
I thought this forum would be an appropriate place for the question and I am really thankful for all the insights shared here.

1 Like

Thank you!!
I stumbled upon this problem when solving for some eigenvalues. For certain choice of parameters, when the values should be same, I was getting two states for this signed zero feature and that was bothering me.
Thanks for your note of caution.

1 Like

@gyrtle , it was mentioned upthread you may be looking for a “one line solution” or some such to resolve this issue. And Julia and Python based program responses were shown as well.

Well, I’m unsure what kind of suggestions you seek but note the dynamic environments such as Julia and Python have minions and minions putting in countless hours pursuing essentially a “monotheistic” path (e.g., have no worries about consensus in a ISO / IEC standard committee where every member can in principle come from a different ideology) in various library solutions.

If you want to match the program behavior of those solutions in Fortran, you can adopt a similar path and “reinvent the wheel” entirely with complex arithmetic and pursue something like the following, for after all, that’s what a language du jour brings, yet another reinvention!

PROGRAM sZ
use complex_m
REAL(DP),PARAMETER :: r = 0.2_dp
type(complex_t(WP=DP)) :: c0
c0 = (0.0,0.5) ; print *, sqrt(c0**2 - r)
c0 = (-0.0,0.5) ; print *, sqrt(c0**2 - r)
END PROGRAM sZ

to get the output you seek:

C:\temp>ifort /standard-semantics sz.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.5.0 Build 20211109_000000

Microsoft (R) Incremental Linker Version 14.30.30706.0

-out:sz.exe
-subsystem:console
sz.obj

C:\temp>sz.exe
(0.00000000000000,0.670820393249937)
(0.00000000000000,0.670820393249937)

C:\temp>

You will note the above is served by this:

Click to see the silly "reinvention of the wheel"
module complex_m
! Add support for other precision to continue the silliness
integer, parameter :: SP = kind(1.0)
integer, parameter :: DP = kind(1D0)
real(DP), parameter :: ZERO = 0.0_dp
real(DP), parameter :: ONE = 1.0_dp
real(DP), parameter :: TWO = 2.0_dp
real(DP), parameter :: EPS = epsilon(ZERO)
type :: complex_t(WP)
integer, kind :: WP = DP
real(kind=WP) :: re = ZERO
real(kind=WP) :: im = ZERO
end type
interface assignment(=)
module procedure assign_complex_dp
end interface
interface operator(**)
module procedure exp_complex_dp
end interface
interface operator(-)
module procedure subtract_complex_real_dp
end interface
interface sqrt
module procedure sqrt_complex_dp
end interface
interface write(formatted)
module procedure write_complex_dp
end interface
private :: assign_complex_dp, exp_complex_dp, subtract_complex_real_dp, sqrt_complex_dp, &
write_complex_dp
contains
subroutine assign_complex_dp( lhs, rhs )
type(complex_t(WP=DP)), intent(out) :: lhs
complex(kind=SP), intent(in)        :: rhs
lhs%re = rhs%re ; lhs%im = rhs%im
end subroutine
function exp_complex_dp( c, n ) result(r)
type(complex_t(WP=DP)), intent(in) :: c
integer, intent(in) :: n
type(complex_t(WP=DP)) :: r
select case ( n )
case ( 2 )
r%re = -( c%re**2 + c%im**2 )
if ( abs(c%re) > EPS ) then
r%im = -TWO*c%re*c%im
else
r%im = ZERO
end if
case default
error stop "exponent not yet supported."
end select
end function
function subtract_complex_real_dp( c, rhs ) result(r)
type(complex_t(WP=DP)), intent(in) :: c
real(kind=DP), intent(in) :: rhs
type(complex_t(WP=DP)) :: r
r%re = c%re - rhs
end function
function sqrt_complex_dp( c ) result(r)
type(complex_t(WP=DP)), intent(in) :: c
type(complex_t(WP=DP)) :: r
if ( abs(c%im) < EPS ) then
r%re = ZERO
r%im = sqrt(abs(c%re))
else
r%re = sqrt( (sqrt(c%re**2+c%im**2) + c%re)/TWO )
r%im = sign(ONE, c%im)*sqrt( (sqrt(c%re**2+c%im**2) - c%re)/TWO )
end if
end function
subroutine write_complex_dp( dtv, lun, iotype, vlist, istat, imsg )

! Argument list
class(complex_t(WP=DP)), intent(in) :: dtv
integer, intent(in)                 :: lun
character(len=*), intent(in)        :: iotype
integer, intent(in)                 :: vlist(:)
integer, intent(out)                :: istat
character (len=*), intent(inout)    :: imsg

! local variable
character(len=20) :: pfmt
integer :: listv
complex(kind=DP) :: c

istat = 0
c%re = dtv%re ; c%im = dtv%im

if ( iotype == "LISTDIRECTED" ) then
write(lun, fmt=*, iostat=istat, iomsg=imsg) c
else if ( iotype(1:2) == "DT" ) then
! vlist(1) is to be used as the field widths of the
! component of the derived type variable. First set up the format to
! be used for output.
if ( size(vlist) > 0 ) then
if ( vlist(1) > 0 ) listv = vlist(1)
end if
write(pfmt,"(*(g0))" ) "(g", listv, ")"
write(lun, fmt=pfmt, advance="no", iostat=istat, iomsg=imsg) c
else
! Not supported
istat = 1
imsg = "Namelist option is not yet supported."
end if

return

end subroutine write_complex_dp
end module
P.S.> Apologies for any errors, it was slapped together rather hastily
1 Like

Well, that ‘silly’ reinvention seems to be quite advanced for me. Thank you for putting it together. I wish I would also be able to slap things together. Someday… 1 Like

I’m unsure what is the problem that needs to be resolved. As I’m sure you know, every nonzero number, real or complex, has two square roots that differ by a sign. You must decide which of those fits your problem. The signed zero issue may be just a distraction, the real question is always which root you want in your application. It is always just a convention, or an implementation detail, which value is returned by the sqrt() function in any language.

Going further, there are always three complex cube roots of a nonzero number, four fourth roots, five fifth roots, and so on. When these arise in your application, you always must choose which value you want.

2 Likes