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:

       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)

The sign of the imaginary part flips.


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


Anton Shterenlikht has a paper on this, 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!

   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)

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
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.30.30706.0
Copyright (C) Microsoft Corporation. All rights reserved.




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, &
   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
               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))
         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
         ! Not supported
         istat = 1
         imsg = "Namelist option is not yet supported."
      end if 
   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… :slight_smile:

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.