Test double complex division in Fortran

test_zcomplexabs.f

*> \brief zabs tests the robustness and precision of the intrinsic ABS for double complex 
*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
*
*> \verbatim
*>
*> Real values for test:
*> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1. Stop on the first success.
*>     Mind that not all platforms might implement subnormal numbers.
*> (2) x = 2**m, where m = MINEXPONENT, ..., 0. Stop on the first success.
*> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV.
*> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1. Stop on the first success.
*>
*> Tests:
*> (a) y = x + 0 * I, |y| = x
*> (b) y = 0 + x * I, |y| = x
*> (c) y = (3/4)*x + x * I, |y| = (5/4)*x whenever (3/4)*x and (5/4)*x can be exactly stored
*> (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2) whenever (1/2)*x can be exactly stored
*>
*> Special cases:
*>
*> (i) Inf propagation
*>    (1) y = Inf + 0 * I, |y| is Inf.
*>    (2) y =-Inf + 0 * I, |y| is Inf.
*>    (3) y = 0 + Inf * I, |y| is Inf.
*>    (4) y = 0 - Inf * I, |y| is Inf.
*>    (5) y = Inf + Inf * I, |y| is Inf.
*>
*> (n) NaN propagation
*>    (1) y = NaN + 0 * I, |y| is NaN.
*>    (2) y = 0 + NaN * I, |y| is NaN.
*>    (3) y = NaN + NaN * I, |y| is NaN.
*>
*> \endverbatim
*
      program zabs

      integer           N, i, nNaN, nInf, min, Max, m
      parameter       ( N = 4, nNaN = 3, nInf = 5 )

      double precision  X( N ), R, threeFourth, fiveFourth, answerC,
     $                  answerD, oneHalf, aInf, aNaN, relDiff, b,
     $                  eps, blueMin, blueMax, Xj, stepX(N), limX(N)
      parameter       ( threeFourth = 3.0d0 / 4,
     $                  fiveFourth = 5.0d0 / 4,
     $                  oneHalf = 1.0d0 / 2 )

      double complex    Y, cInf( nInf ), cNaN( nNaN )
      intrinsic         ABS, DBLE, RADIX, CEILING, TINY, DIGITS, SQRT,
     $                  MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX,
     $                  EPSILON
*
      min = MINEXPONENT(0.0d0)
      Max = MAXEXPONENT(0.0d0)
      m = DIGITS(0.0d0)
      b = DBLE(RADIX(0.0d0))
      eps = EPSILON(0.0d0)
      blueMin = b**CEILING( (min - 1) * 0.5d0 )
      blueMax = b**FLOOR( (Max - m + 1) * 0.5d0 )
*
      X(1) = TINY(0.0d0) * b**( DBLE(1-m) )
      X(2) = TINY(0.0d0)
      X(3) = HUGE(0.0d0)
      X(4) = b**( DBLE(Max-1) )
*
      stepX(1) = 2.0
      stepX(2) = 2.0
      stepX(3) = 0.0
      stepX(4) = 0.5
*
      limX(1) = X(2)
      limX(2) = 1.0
      limX(3) = 0.0
      limX(4) = 2.0
*
      print *, '# X :=', X
      print *, '# Blue min constant :=', blueMin
      print *, '# Blue max constant :=', blueMax
*
      Xj = X(1)
      if( Xj .eq. 0.0d0 ) then
        print *, "# Subnormal numbers treated as 0"
      else
        do 100 i = 1, N
            Xj = X(i)
            if( Xj .eq. 0.0d0 ) print *, 
     $         "# Subnormal numbers may be treated as 0"
 100    continue
      endif
*
      aInf = X(3) * 2
      cInf(1) = DCMPLX( aInf, 0.0d0 )
      cInf(2) = DCMPLX(-aInf, 0.0d0 )
      cInf(3) = DCMPLX( 0.0d0, aInf )
      cInf(4) = DCMPLX( 0.0d0,-aInf )
      cInf(5) = DCMPLX( aInf,  aInf )
*
      aNaN = aInf / aInf
      cNaN(1) = DCMPLX( aNaN, 0.0d0 )
      cNaN(2) = DCMPLX( 0.0d0, aNaN )
      cNaN(3) = DCMPLX( aNaN,  aNaN )
*
*     Test (a) y = x + 0 * I, |y| = x
      do 10 i = 1, N
        Xj = X(i)
        if( Xj .eq. 0.0d0 ) then
            print *, "# [a] Subnormal numbers may be treated as 0"
        else
            do while( Xj .ne. limX(i) )
                  Y = DCMPLX( Xj, 0.0d0 )
                  R = ABS( Y )
                  if( R .ne. Xj ) then
                        relDiff = ABS(R-Xj)/Xj
                        WRITE( *, FMT = 9999 ) 'a',i, Y, R, Xj, relDiff
                  endif
                  Xj = Xj * stepX(i)
            end do
        endif
  10  continue
*
*     Test (b) y = 0 + x * I, |y| = x
      do 20 i = 1, N
        Xj = X(i)
        if( Xj .eq. 0.0d0 ) then
            print *, "# [b] Subnormal numbers may be treated as 0"
        else
            do while( Xj .ne. limX(i) )
                  Y = DCMPLX( 0.0d0, Xj )
                  R = ABS( Y )
                  if( R .ne. Xj ) then
                        relDiff = ABS(R-Xj)/Xj
                        WRITE( *, FMT = 9999 ) 'b',i, Y, R, Xj, relDiff
                  endif
                  Xj = Xj * stepX(i)
            end do
        endif
  20  continue
*
*     Test (c) y = (3/4)*x + x * I, |y| = (5/4)*x
      do 30 i = 1, N
        if( i .eq. 3 ) go to 30
        if( i .eq. 1 ) then
            Xj = 4*X(i)
        else
            Xj = X(i)
        endif
        if( Xj .eq. 0.0d0 ) then
            print *, "# [c] Subnormal numbers may be treated as 0"
        else
            do while( Xj .ne. limX(i) )
                  answerC = fiveFourth * Xj
                  Y = DCMPLX( threeFourth * Xj, Xj )
                  R = ABS( Y )
                  if( R .ne. answerC ) then
                        relDiff = ABS(R-answerC)/answerC
                        WRITE( *, FMT = 9999 ) 'c',i, Y, R, answerC, relDiff
                  endif
                  Xj = Xj * stepX(i)
            end do
        endif
  30  continue
*
*     Test (d) y = (1/2)*x + (1/2)*x * I, |y| = (1/2)*x*sqrt(2)
      do 40 i = 1, N
        if( i .eq. 1 ) then
            Xj = 2*X(i)
        else
            Xj = X(i)
        endif
        if( Xj .eq. 0.0d0 ) then
            print *, "# [d] Subnormal numbers may be treated as 0"
        else
            do while( Xj .ne. limX(i) )
                  answerD = (oneHalf * Xj) * SQRT(2.0d0)
                  if( answerD .eq. 0.0d0 ) then
                  print *, "# [d] Subnormal numbers may be treated as 0"
                  else
                  Y = DCMPLX( oneHalf * Xj, oneHalf * Xj )
                  R = ABS( Y )
                  relDiff = ABS(R-answerD)/answerD
                  if( relDiff .ge. (0.5*eps) )
     $                  WRITE( *, FMT = 9999 ) 'd',i, Y, R, answerD, relDiff
                  endif
                  Xj = Xj * stepX(i)
            end do
        endif 
  40  continue
*
*     Test (e) Infs
      do 50 i = 1, nInf
        Y = cInf(i)
        R = ABS( Y )
        if( .not.(R .gt. HUGE(0.0d0)) ) then
            WRITE( *, FMT = 9999 ) 'i',i, Y, R, aInf, aInf
        endif
  50  continue
*
*     Test (f) NaNs
      do 60 i = 1, nNaN
        Y = cNaN(i)
        R = ABS( Y )
        if( R .eq. R ) then
            WRITE( *, FMT = 9998 ) 'n',i, Y, R
        endif
  60  continue
*
 9998 FORMAT( '[',A1,I1, '] ABS( ', (ES10.3,SP,ES10.3,"*I"), ' ) = ',
     $        ES10.3, ' differs from NaN' )
*
 9999 FORMAT( '[',A1,I1, '] ABS( ', (ES10.3,SP,ES10.3,"*I"), ' ) = ',
     $  ES10.3,' differs from ', ES10.3, ' #relative diff = ', ES10.3 )
*
*     End of zabs
*
      END