Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28386 Discussions

Possible bug in exponentiation complex type

Kyle_T_
Beginner
369 Views

I believe I may have discovered a bug in the ifort version 15.0.3 for Mac OSX regarding exponentiating a number of type complex.  The following code evaluates

(a + bi)^(1 + 0i):

program test

  implicit none

  integer, parameter :: dp=selected_real_kind(15, 307)

  complex(dp) :: quantity
  complex(dp) :: power
  complex(dp) :: computed_result

  real(dp), parameter :: zero = 0.0_dp
  real(dp), parameter :: half = 0.5_dp
  real(dp), parameter :: one  = 1.0_dp

continue

  ! Values that produce correct result
  !quantity = (0.5_dp,0.5_dp)
  !power    = (1._dp,0._dp)

  ! Values that produce the incorrect result
  quantity = (-2.183392967490150E-003_dp,-2.183392967490262E-019_dp)
  power    = (one,zero)

  computed_result = (quantity)**power

  write(*,*) "original value = ",quantity
  write(*,*) "computed value = ",computed_result

  if (aimag(quantity) /= aimag(computed_result)) then
    write(*,*) "FAIL: original value and computed value are not equal"
  else
    write(*,*) "PASS: original value and computed value are equal"
  end if

end program test

I would expect that exponentiating a complex number by one would always return the original number, but that does not appear to be the case here.  The issue remains regardless of the precision used, but does appear to be dependent on the magnitude of the complex quantity.  Is this a bug in the compiler, or am I misinterpreting how the expression is evaluated?

0 Kudos
1 Solution
jimdempseyatthecove
Honored Contributor III
369 Views

Your numbers for quantity have a magnitude difference of E-16. The precision of DP is ~15.95 digits. Your expression will therefore present rounding errors (to the point of total truncation). You may have to special case the power=(one,zero) condition if this is a requirement for your use.

Jim Dempsey

View solution in original post

0 Kudos
4 Replies
jimdempseyatthecove
Honored Contributor III
370 Views

Your numbers for quantity have a magnitude difference of E-16. The precision of DP is ~15.95 digits. Your expression will therefore present rounding errors (to the point of total truncation). You may have to special case the power=(one,zero) condition if this is a requirement for your use.

Jim Dempsey

0 Kudos
TimP
Honored Contributor III
369 Views

If you want full accuracy for  **1 you should use the integer exponentiation.  Fortran standard used to permit **0 to be undefined; I don't know if that is still so.  It's a reasonable expectation for it to work.

ifort often treats constant real integral exponents as integer exponents, but that can lead to failure to detect non-standard code.  There probably is no short-cutting for complex exponents with 0 imaginary part, nor any special treatment to protect results with extra precision such as one would hope for with real scalar exponentiation (by splitting the exponent into integral and fractional parts).  So it's not surprising if the result is no better than you would get with a naive implementation.

The classic celefunt test suite might interest you.

0 Kudos
Kyle_T_
Beginner
369 Views

Thank you both for the prompt response. Jim, you are exactly correct.  The issue is indeed truncation error, and I can confirm that reducing the magnitude difference fixes (mitigates) the issue.

Tim, excellent suggestion on using integer exponentation for the **1 case.  That does indeed recover the full precision, and it may be suitable for my needs.  Also, thank you for pointing me to the classic celefunt test suite.  I found it quite interesting and useful.

0 Kudos
jimdempseyatthecove
Honored Contributor III
369 Views

Kyle,

An easy way to test for integer values is to use IBITS to extract the mantissa of the given floating point value. If the mantissa (with implicit 1 bit) is 0, then the IEEE-754 encoded floating point number is an integer (or 0),It could also be one of the NaN's. That would be the quickest way to make a test: IF(testForInteger(x)) THEN ... ELSE ... ENDIF

Jim Dempsey

0 Kudos
Reply