Community
cancel
Showing results for
Did you mean:
Highlighted
Beginner

## Possible bug in exponentiation complex type

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?

Tags (1)

Accepted Solutions
Highlighted
Black Belt

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

4 Replies
Highlighted
Black Belt
1 View

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

Highlighted
Black Belt

## If you want full accuracy for

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.

Highlighted
Beginner

## Thank you both for the prompt

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.

Highlighted
Black Belt