!     COMPLEX.F90  Copyright(c) 1994 - 1998, Lahey Computer Systems, Inc.
!     Copying for sale requires permission from Lahey Computers Systems.
!     Otherwise, distribution of all or part of this file is permitted
!     if these four lines are included.

!  Complex derived-type MODULE with some of the overloaded operators
   MODULE complex
     IMPLICIT NONE            !  All variables' type must be declared
     PRIVATE                  !  Only PUBLIC statement names global

!  PUBLIC, make names and operators available globally (to USErs of the MODULE)
     PUBLIC  ::  cx, fpkind, ASSIGNMENT(=), OPERATOR(+), OPERATOR(-),        &
                             OPERATOR(*),   OPERATOR(/), OPERATOR(**)

     INTEGER, PARAMETER  ::  fpkind = KIND(0.0D0)

     TYPE cx                  !  Can't use type names (REAL, COMPLEX, ...)
       REAL(fpkind) :: real, imag
     END TYPE cx

!  These 6 INTERFACE definitions overload the intrinsic binary operators

!  Code for the unary operators is not written, needed for completeness

     INTERFACE ASSIGNMENT  (=)
       MODULE PROCEDURE cx_c_asg_i                 !  Need:  cx_c_asg_r
     END INTERFACE

     INTERFACE OPERATOR  (+)
       MODULE PROCEDURE cx_cadr, cx_radc, cx_cadc  !  Need:  cx_cadi, cx_iadc
     END INTERFACE

     INTERFACE OPERATOR  (-)
       MODULE PROCEDURE cx_csbc    !  Need: cx_csbi, cx_isbc, cx_csbr, cx_rsbc
     END INTERFACE

     INTERFACE OPERATOR  (*)
       MODULE PROCEDURE cx_cmlc    !  Need: cx_cmli, cx_imlc, cx_cmlr, cx_rmlc
     END INTERFACE

     INTERFACE OPERATOR  (/)
       MODULE PROCEDURE cx_cdvc    !  Need: cx_idvc, cx_cdvi, cx_cdvr, cx_rdvc
     END INTERFACE

     INTERFACE OPERATOR  (**)
       MODULE PROCEDURE cx_cxpi    !  Need: cx_ixpc, cx_cxpr, cx_rxpc, cx_cxpc
     END INTERFACE

     CONTAINS

!Assignments.  Missing cx = REALs
     SUBROUTINE cx_c_asg_i(z, i)   !  complex_target = integer_value
       INTEGER, INTENT(IN) :: i;        TYPE (cx), INTENT(OUT) :: z

       z%real = i;                      z%imag = 0.0
     END SUBROUTINE cx_c_asg_i

!Addition
     FUNCTION cx_cadr(z, r)        !  complex + real
     ENTRY cx_radc(r, z)           !  real + complex
       TYPE (cx)             :: cx_cadr, cx_radc
       TYPE (cx), INTENT(IN) :: z;    REAL(fpkind), INTENT(IN) :: r

       cx_cadr%real = z%real + r;    cx_cadr%imag = z%imag
     END FUNCTION cx_cadr
 
     FUNCTION cx_cadc(zl, zr)      !  complex + complex
       TYPE (cx) :: cx_cadc;            TYPE (cx), INTENT(IN) :: zl, zr

       cx_cadc%real = zl%real + zr%real
       cx_cadc%imag = zl%imag + zr%imag
     END FUNCTION cx_cadc


!Subtraction
     FUNCTION cx_csbc(zl, zr)      !  complex - complex
       TYPE (cx)             :: cx_csbc;   TYPE (cx), INTENT(IN) ::  zl, zr

       cx_csbc%real = zl%real - zr%real
       cx_csbc%imag = zl%imag - zr%imag
     END FUNCTION cx_csbc


!Multiplication
     FUNCTION cx_cmlc(zl, zr)      !  complex*complex
       TYPE (cx)             :: cx_cmlc;   TYPE (cx), INTENT(IN) ::  zl, zr

       cx_cmlc%real = zl%real*zr%real - zl%imag*zr%imag
       cx_cmlc%imag = zl%real*zr%imag + zl%imag*zr%real
     END FUNCTION cx_cmlc


!Division
     FUNCTION cx_cdvc(zl, zr)      !  complex_left/complex_right
       REAL (KIND = KIND (0.0D0)) :: temp
       TYPE (cx) :: cx_cdvc;    TYPE (cx), INTENT(IN) :: zl, zr

       IF ( zr%real == 0.0 ) THEN
!  Divisor, zr, real part is 0.0, imaginary part could be 0.0
         IF ( zr%imag == 0.0 )  STOP                  &
            "ABORT:  Complex division, divisor == (0.0, 0.0)"
!  Divisor, zr, real part is 0.0, imaginary part is nonzero
         cx_cdvc%real = zl%imag/zr%imag
         cx_cdvc%imag = zl%real*zr%imag

       ELSE IF ( zr%imag == 0.0 )  THEN
!  Divisor,zr, real part is nonzero, imaginary part is 0.0
         cx_cdvc%real = zl%real/zr%real
         cx_cdvc%imag = zl%imag*zr%real

       ELSE                 !  Both parts of the divisor, zr, are nonzero
         temp = zr%real**2D0 + zr%imag**2D0!     DP divisor, squares sum
         cx_cdvc%real = (zl%real*zr%real + zl%imag*zr%imag)/temp
         cx_cdvc%imag = (zl%imag*zr%real - zl%real*zr%imag)/temp
       END IF
     END FUNCTION cx_cdvc
 
!Exponentiation
     FUNCTION cx_cxpi(z, i)              !  complex**integer
       INTEGER, INTENT(IN)  ::  i
       INTEGER              ::  in, n
       TYPE (cx)            ::  cx_cxpi, z_xp2_n;   TYPE (cx), INTENT(IN) :: z

!  Special case when z part(s)is zero; ABORT or faster execution
       IF ( z%real == 0.0  .AND.  z%imag == 0.0 )  THEN
         IF ( i == 0 )  STOP "ABORT:  cx(0.0, 0.0)**0 is not defined!"
         IF ( i < 0 )   STOP "ABORT:  cx(0.0, 0.0) to a negative power"
         cx_cxpi = cx(0.0, 0.0)        !  zero ** nonzero is zero
         RETURN
       END IF

       in = ABS(i)                     !  Initialize local variable

       IF ( IAND(in, 1) /= 0 ) THEN
         cx_cxpi = z                   !  n odd.  Fortran can do =
       ELSE
         cx_cxpi = cx(1.0, 0.0)        !  n even.  non-zero**0 is 1.0
       END IF

       n = 2                           !  Loop initialization
       z_xp2_n = z                     !  Fortran can do =

!  Loop forms the product of the powers of two that sum to ABS(i)
       DO WHILE ( in >= n ) 
         z_xp2_n = z_xp2_n*z_xp2_n                         !  * Overload
         IF (IAND(in, n) /= 0)  cx_cxpi = cx_cxpi*z_xp2_n  !  * Overload
         n = 2*n
       END DO

!  For i negative, reciprocate
     IF ( i < 0 )  cx_cxpi = cx(1.0, 0.0)/cx_cxpi          !  / Overload
     END FUNCTION cx_cxpi

   END MODULE complex