1/z when z is a complex number with a +Inf component

1/z when z is a complex number with a +Inf component

Dear Fortran experts, As a sub problem in my code I need to evaluate the following expression: lim 1 x -> 0+ ------------------ sqrt(-1) + 1/x The expression above is well defined in the extended complex plane. 1/x ---> Inf as x ---> 0. sqrt(-1) + Inf = Inf. 1/Inf = 0. To try it out I wrote a code, compiled and ran it with three different Fortran 95 processors and got two different results. REAL X COMPLEX I PARAMETER(I=(0.0,1.0)) INTRINSIC SIGN X=SIGN(0.0,+1.0) PRINT*,'X=',X PRINT*,'1/X=',1/X PRINT*,'I+1/X=',I+1/X PRINT*,'1/(I+1/X)=',1/(I+1/X) X=SIGN(0.0,-1.0) PRINT*,'X=',X PRINT*,'1/X=',1/X PRINT*,'I+1/X=',I+1/X PRINT*,'1/(I+1/X)=',1/(I+1/X) END In GFortran and Open64 (I+1/X) evaluates to (0.0,1.0) + (+Inf,0.0) = (+Inf,1.0) (assuming X=+0.0) which [probably?] should be interpreted as a the complex infinity. Intel Fortran I got (NaN,0.0) as result. Is this a bug in Intel Fortran or is there something here which I am not thinking of? Am I violating a Fortran rule or IEEE arithmetic rule somewhere? Dr. Fortran (Steve Lionel) kindly pointed it out to me that IEEE 754 does not define complex arithmetic. I assume then that Fortran has at least *some* rules about how complex arithmetic on the *extended* complex plane should be executed by the processor? Keywords: Fortran 95, IEEE 754, complex, extended real, extended complex, arithmetic, reciprocal.

19 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

Quote:

Henrik H. wrote:

[snip]

Intel Fortran I got (NaN,0.0) as result.

A typo there. The (NaN,0.0) result is for the *entire* expression.

I better attach the output so you can see yourself what is going on. I am terrible bad at explaing things. ;-)

Intel output:

Quote:

 X=  0.0000000E+00
 1/X= Infinity
 I+1/X= (Infinity,1.000000)
 1/(I+1/X)= (NaN,0.0000000E+00)
 X=  0.0000000E+00
 1/X= -Infinity
 I+1/X= (-Infinity,1.000000)
 1/(I+1/X)= (NaN,0.0000000E+00)

 

GFortran output:

Quote:

 X=   0.0000000
 1/X=       +Infinity
 I+1/X= (      +Infinity,  1.0000000    )
 1/(I+1/X)= (  0.0000000    ,  0.0000000    )
 X=  -0.0000000
 1/X=       -Infinity
 I+1/X= (      -Infinity,  1.0000000    )
 1/(I+1/X)= ( -0.0000000    , -0.0000000    )

In order to have treatment of division in ifort similar to what gfortran gives when you don't set -ffast-math or -freciprocal-math, you should set corresponding ifort options -prec-div -prec-sqrt (or options such as -fp-model source which imply those).  

I'm not certain I could extract your test case correctly from what you posted.  It looks like it may involve compile-time evaluation, which current gfortran carries out by gnu multiple precision math library.  This doesn't necessarily produce the same result you would get by run-time evaluation.

If your question is associated with the more aggressive default options of commercial compilers, opposed to what gfortran does as the result of consensus from users, it's certainly a more difficult question than it should be, one which will occupy some space if I am able to publish on it.

Believe it or not, the Fortran standard doesn't even say 2+2=4. It does not require use of IEEE 754 arithmetic and doesn't establish any interpretation other than the generic mathematical one. For example, the sum total definition of the / operator is "Divide x1 by x2".

That said, I can't reproduce exactly the results you claim for Intel Fortran. I will note that you want to use the "-assume minus0" option to get the behavior of -0 you seem to want. When I use this, I get:

 X=  0.0000000E+00
 1/X=       Infinity
 I+1/X= (Infinity,1.000000)
 1/(I+1/X)= (NaN,NaN)
 X= -0.0000000E+00
 1/X=      -Infinity
 I+1/X= (-Infinity,1.000000)
 1/(I+1/X)= (NaN,NaN)

I was going to suggest also using -fp-model strict, but in my experiments that didn't change the results. Not a bad idea, though, if you are being picky. The only divergence seems to be the last expression. Complex divide is a complex (!) operation - especially when infinities are involved. That said, the NaN result seems a bit off to me and I'll ask one of my math library colleagues for an opinion.

Steve - Intel Developer Support

Ok, I've studied this a bit more and coded some test cases. I believe that Intel Fortran is correct here.

I looked at the definition of complex divide at Wikipedia. Let's take the first case with positive infinity.

Numerator = (1,0)
Denominator = (Inf,1)

Going by the notation in the article, the numerator is a+bi and the denominator is c+di. This gives us:

a = 1
b = 0
c = Inf
d = 1

The real part of the quotient is computed as ((a*c)+(b*d))/((c*c)+(d*d)) So let's see what we get:

a*c is 1*Inf or Inf
b*d is 0*1 or 0
c*c is Inf*Inf or Inf
d*d is 1*1 or 1

So we now have (Inf+0)/(Inf+1) or Inf/Inf. That's a NaN for the real part

The imaginary part is computed as ((b*c)-(a*d))/((c*c)+(d*d)) 

b*c is 0*Inf or NaN
a*d is 1*1 or 1

So for the imaginary part we have (NaN-1)/(Inf+1) or NaN/Inf, which gives us NaN for the imaginary part.

I could have made an error here - if so, let me know.

Steve - Intel Developer Support

According to http://en.wikipedia.org/wiki/Riemann_sphere

    z/Inf = 0

for any finite complex number z. Why does the formula Steve posted above not work? Lets see how it's usually done in the text books, for finite z (|z|<Inf):

Let z = a + b i where i = sqrt(-1). It hold that for any extended complex number z:

   1/z = 1/(a+b i) = (1/(a+b i)) * 1

The next step is, usually, that you rewrite the rightmost factor 1 as 1 = (a-b i)/(a - b i) the complex conjugate of z. However, I think that this is an invalid operation. Since z = Inf in the extended complex plane we are trying to use 1 = Inf/Inf which of course is equal to NaN, not 1.

If this is the correct way to interpret it from the Fortran standard I don't know. But from a mathematical perspective it must be right to evaluate 1/z = 0 in the case when z = Inf (note that there is only one Infinity in the extended complex plane).

We could also consider the limit value of 1/z as z ---> Inf: By expressing z in polar coordinates as z = r e^(i theta) we get that 1/z is (absolutely) convergent with limit 0 outside a arbitrary large disc { z : |z| > R}:

   |1/z| = |(1/r) exp(-i theta)| = 1/r |exp(-i theta| < 1/R ---> 0 as R ---> Inf,

or with the formula posted in the previous post:

   |1/z| = |(1/z) * conj(z)/conj(z)| = |conj(z)| / |z|^2 = (r exp(-i theta)) / r^2 <= (1/r) |exp(-i theta)| < 1/R ---> 0 as R ---> Inf.

I was thinking of a reasonable way to determine if z=Inf (complex infinity) or z=NaN (complex NaN) such that it does not matter if w=(Inf,1.0), w=(2.0,-Inf), ... or w=(Inf,Inf) they are all equivalent to the extended complex Infinity and the same for the complex NaN.

A function computing reciprocal values of extended complex numbers could be of the form:

[code]

complex elemental function reciprocal(z)
  complex, intent(in) :: z
  if (isNaN(z)) then
    reciprocal = z
    return
  endif
  if (isInf(z)) then
    reciprocal=(0.0,0.0) !! I don't know if we should or even can care about signed zeros here
    return
  endif
  reciprocal=1/z !! Now we trust the Fortran processor to do the reciprocal computation itself.
contains
  logical elemental function isNaN(z)
    complex, intent(in) :: z
    isNaN = real(z) /= real(z) .or. aimag(z) /= aimag(z)
  end function isNaN
  logical elemental function isInf(z)
    complex, intent(in) :: z
    real, parameter :: realmax = huge(1.0)
    isInf = abs(real(z)) > realmax .or. abs(aimag(z)) > realmax
  end function isInf
end function reciprocal

[/code]

Again, the Fortran standard is entirely silent on this. But IEEE 754 defines division by infinity only if the divisor is finite, otherwise it's "indeterminate". I am not a mathematician and would fall off a Reimann sphere, but every reference I could find on arithmetic with infinities had inf/inf not being finite.

Steve - Intel Developer Support

Quote:

Tim Prince wrote:

In order to have treatment of division in ifort similar to what gfortran gives when you don't set -ffast-math or -freciprocal-math, you should set corresponding ifort options -prec-div -prec-sqrt (or options such as -fp-model source which imply those).  

I'm not certain I could extract your test case correctly from what you posted.  It looks like it may involve compile-time evaluation, which current gfortran carries out by gnu multiple precision math library.  This doesn't necessarily produce the same result you would get by run-time evaluation.

If your question is associated with the more aggressive default options of commercial compilers, opposed to what gfortran does as the result of consensus from users, it's certainly a more difficult question than it should be, one which will occupy some space if I am able to publish on it.

The option "/fp:strict" did change the end result. However I am still not convinced it is a "expected" result.

Quote:

$ ./expr
 X=  0.0000000E+00
 1/X= Infinity
 I+1/X= (Infinity,1.000000)
 1/(I+1/X)= (NaN,NaN)
 X=  0.0000000E+00
 1/X= -Infinity
 I+1/X= (-Infinity,1.000000)
 1/(I+1/X)= (NaN,NaN)

If we compare this to a REAL computation of the same kind we get 0.0. As I expected.

Quote:

$ cat expr2.f
      REAL X
      INTRINSIC SIGN
      X=SIGN(0.0,+1.0)
      PRINT*,'X=',X
      PRINT*,'1/X=',1/X
      PRINT*,'1+1/X=',1+1/X
      PRINT*,'1/(1+1/X)=',1/(1+1/X)
      X=SIGN(0.0,-1.0)
      PRINT*,'X=',X
      PRINT*,'1/X=',1/X
      PRINT*,'1+1/X=',1+1/X
      PRINT*,'1/(1+1/X)=',1/(1+1/X)
      END

$ ifort /fp:strict expr2.f
Intel(R) Visual Fortran Intel(R) 64 Compiler XE for applications running on Intel(R) 64, Version 12.0.2.154 Build 20110112
Copyright (C) 1985-2011 Intel Corporation.  All rights reserved.

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

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

$ ./expr2
 X=  0.0000000E+00
 1/X= Infinity
 1+1/X= Infinity
 1/(1+1/X)=  0.0000000E+00
 X=  0.0000000E+00
 1/X= -Infinity
 1+1/X= -Infinity
 1/(1+1/X)=  0.0000000E+00

The real version is not at all the same computation. 1/Inf is indeed 0. But Inf/Inf is NaN.

Here's how I modeled the complex computation.

    program Console2

	    

	    use, intrinsic :: ieee_arithmetic
    implicit none
    complex numerator, divisor

	    real a,b,c,d

	    real ac, bd, bc, ad, c2,d2

	    

	    a = 1.

	    b = 0.

	    c = IEEE_VALUE(0.0,IEEE_POSITIVE_INF)

	    d = 1

	    

	    

	    numerator = cmplx(a,b)

	    divisor = cmplx(c,d)

	    

	    print *, "a, b, c, d: ",a,b,c,d

	    print *, "numerator, divisor: ", numerator, divisor

	    

	    ac = a*c

	    bd = b*d

	    bc = b*c

	    ad = a*d

	    c2 = c*c

	    d2 = d*d

	    

	    print *, "ac, bd, bc, ad: ", ac,bd,bc,ad

	    print *, "c**2, d**2: ", c2, d2

	    

	    print *, "ac+bd=", ac+bd

	    print *, "c**2+d**2", c2+d2

	    print *, "Quotient (real part):", (ac+bd)/(c2+d2)

	    

	    print *, "bc-ad=", bc-ad

	    print *, "c**2+d**2", c2+d2

	    print *, "Quotient (imaginary part):", (bc-ad)/(c2+d2)

	    end program Console2

	

C:\Projects>ifort /fp:strict console2.f90
Intel(R) Visual Fortran Compiler XE for applications running on IA-32, Version 14.0.1.139 Build 20131008
Copyright (C) 1985-2013 Intel Corporation.  All rights reserved.

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

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

C:\Projects>console2.exe
 a, b, c, d:    1.000000      0.0000000E+00       Infinity   1.000000
 numerator, divisor:  (1.000000,0.0000000E+00) (Infinity,1.000000)
 ac, bd, bc, ad:        Infinity  0.0000000E+00            NaN   1.000000
 c**2, d**2:        Infinity   1.000000
 ac+bd=       Infinity
 c**2+d**2       Infinity
 Quotient (real part):            NaN
 bc-ad=            NaN
 c**2+d**2       Infinity
 Quotient (imaginary part):            NaN

 

Steve - Intel Developer Support

Quote:

Steve Lionel (Intel) wrote:
The real version is not at all the same computation. 1/Inf is indeed 0. But Inf/Inf is NaN.
I appreciate the interesting discussion we are having. Thank you Steve and Tim for your time and interest. Steve: I agreed with both your statements. But I fail to see how the computation done by console2.f90 is relevant. Applying the formula z/w = (z*conj(w))/|w|^2 is only valid when w is a finite complex number. The extended complex plane has (one) infinity, denoted here cInfinity, which satisfy z/cInfinity = 0 for any finite complex number z, including z=1. I think the discussion boils down to this: How should the complex infinity be handled in Fortran? I guess this is uncharted territory. The Fortran standard seems to be (for me unexpectedly) vague on the mathematical correctness of complex arithmetic including complex infinity. Maybe this is nothing else than an oversight in the standard. Or it is intentional. But I cannot see why it would be so. In my opinion, not knowing the result of 1/cInfinity makes Fortran unreliable for computations in (native) complex arithmetic. If we accept "my" definition, isinf(z) := isinf(real(z)).or.isinf(aimag(z)) which is a reasonable mathematically consistent definition, then we should have 1/(Inf,1) == 1/cInfinity which, mathematically and most indisputably, is equal to 0. Maybe we must just agreed to disagree about the subject matter. ;-) I think it would be a bit ashame to do so, because it would lead to a processor dependent behavior which I think the Fortran community gains no benefit from. I wish I still had access to some more Fortran processors I could try to experiment with. If I knew that this was a "known issue" among many processors I would not be so concerned.

As I wrote earlier, the Fortran standard is entirely silent on such matters, deliberately so. There is the IEEE_ARITHMETIC intrinsic module which reflects IEEE 754, but the main part of the standard is agnostic regarding the floating point environment - necessarily so since there are Fortran implementations on processors that don't use IEEE-754 FP. Also, as I wrote earlier, IEEE 754 does not define complex operations.

My take on this is that the mathematical community is generally quite satisfied with the Fortran language, but in the area of complex arithmetic, there are sometimes multiple interpretations or values and an implementation has to pick one.

Steve - Intel Developer Support

Steve, I think your arithmetic is incorrect because it doesn't take into account the possibility that a**2+b**2 may overflow when the final result should not.

program p

	   implicit none

	   integer, parameter :: dp = kind(1.0d0)

	   complex(dp) c

	   real(dp) a, b(2)
   a = sqrt(huge(a))

	   c = (100,100)*a

	   write(*,*) 'c = ',c

	   a = abs(c)

	   write(*,*) 'abs(c) = ', a

	   c = 1/c

	   write(*,*) '1/c = ', c

	   a = sqrt(huge(a))

	   b = [100,100]*a

	   write(*,*) 'b = ',b

	   a = norm2(b)

	   write(*,*) 'norm2(b) = ', norm2(b)

	end program p

Quote:

 c =  (  1.3407807929942595E+156,  1.3407807929942595E+156)
 abs(c) =    1.8961503816218350E+156
 1/c =  (  3.7291703656001041E-157, -3.7291703656001041E-157)
 b =    1.3407807929942595E+156   1.3407807929942595E+156
 norm2(b) =    1.8961503816218350E+156

With your arithmetic above, we would get infinities or zeros rather than the reasonable results displayed.

Default options should protect against intermediate over/underflow in operations such as complex sqrt, abs, divide.  The option -complex-limited-range (implied by -fp-model fast=2) removes that protection.  If RO is making the point that writing these operations out with reals won't necessarily produce the same result, I agree.

Steve's suggested program is a starting point for the reader to use and expand upon. Fortran does not provide for infinite precision, nor does it discriminately promote numeric expressions to larger internal precisions to postpone potential overflow conditions. This is left to the programmer. Note, promoting REAL(4) to REAL(16) would still not be sufficient to keep outlier cases from overflowing.

REAL(4) - Single Precision has 8 bits for exponent, 24 bits of mantissa (23 stored, 1 implicit)
REAL(8) - Double Precision has 11 bits for exponent, 53 bits of mantissa (52 stored, 1 implicit)
REAL(16) - Quadruple Precision has 15 bits for exponent, 113 bits of mantissa (***)
(*** not sure about the implicit 1 bit, would have to research, there is 112 bits of storage available for 112 actual bits of mantissa without implicit 1 or 113 bits of mantissa with implicit 1)

REAL(4) would require a promotion to an internal format having 2x8-bits for exponent, 2x24-bits for mantissa to assure non-overflow
While REAL(8) can avoid (additional) loss in the mantissa, it cannot assure all cases of exponent overflow, neither could promotion to REAL(16). Protecting REAL(4) from all cases of overflow for this implementation of complex division is problematic using provided precisions. Protecting REAL(8) is even more so problematic. Therefore you have to live within your means. The Fortran standards committee is usually silent (publicly silent, though internally I suspect near fisticuffs) on these issues (e.g. should the outlier case result in QNaN, SNaN, Inf, 1, 0).

There arbitrary precision libraries (typically C++) that you could use for the cases where the internal temporary products of non-infinity produce infinity (or near huge).

Jim Dempsey

www.quickthreadprogramming.com

-no-complex-limited-range is carried out for real(4) by promotion of in-line expanded intrinsics to real(8), usually retaining simd optimization, but accepting the reduction in simd lanes (just one complex*16 at a time, so not necessarily considered vectorized for simd-128).  In real(8), it's typically in-line expansion with multiple case branches depending on relative magnitude of real and imaginary, losing vectorization for the expansion of those intrinsics.  -opt-report doesn't show these choices.

linux distributions (and cygwin et al.) nowadays include the multiple precision packages gmp, mpfr, mpc, which, as Jim said, aren't designed for calling directly from Fortran.   I don't know if anyone investigated iso_c_binding for this.  There are some well-known traditional Fortran multiple precision libraries.

gfortran uses those multiple precision packages for compile-time expression evaluation, thus the point I made earlier that the result could be different between compile and run-time expression evaluation, or between different compilers.

Actually, my program was just intended to be a demonstration of the canonical algorithm for complex divide, showing that if you follow the steps you end up with the same result as ifort does for the divide. It is certainly not intended as a general purpose routine

Steve - Intel Developer Support

Steve already mentioned that the Fortran standard does not address the issue discussed here. So it comes down to a "quality of implementation" issue.

On the other hand, the C99 standard specifies complex arithmetic in Annex G.  Binary operators are treated in G.5. Not everything in the specification can be said to be intuitive, and I do not know whether everything is always consistent (e.g. what is being considered "infinite"). However, you may need to use a #pragma CX_LIMITED_RANGE set to "off" to get the specified behavior.

G.5.1, paragraph 4:

The * and / operators satisfy the following infinity properties for all real, imaginary, and
complex operands:315)
— if one operand is an infinity and the other operand is a nonzero finite number or an
   infinity, then the result of the * operator is an infinity;
— if the first operand is an infinity and the second operand is a finite number, then the
   result of the / operator is an infinity;
— if the first operand is a finite number and the second operand is an infinity, then the
   result of the / operator is a zero;
— if the first operand is a nonzero finite number or an infinity and the second operand is
   a zero, then the result of the / operator is an infinity.

I find it somewhat disturbing that in complex arithmetic infinities more or less override NaNs, but it's the C spec. after all.  The question is, does or should Intel Fortran's -no-complex-limited-range behave the same as C99 Annex G?

 

I disagree that this is "quality of implementation". C has some peculiar definitions of arithmetic, some of which explicitly conflict with the Fortran standard. There's nothing to suggest that a Fortran implementation "should" match C99's behavior in any arithmetic operation.

C99 specifies behavior here, though as you say it is a bit odd.  A C99 compiler needs to do what it says. A Fortran implementation does not, and there's nothing that says that using a particular ifort option will necessarily get you C99 semantics.

Steve - Intel Developer Support

Leave a Comment

Please sign in to add a comment. Not a member? Join today