############################################################################
#
#	File:     complex.icn
#
#	Subject:  Procedures to perform complex arithmetic
#
#	Author:   Ralph E. Griswold
#
#	Date:     May 26, 2010
#
############################################################################
#
#   This file is in the public domain.
#
############################################################################
#  
#  The following procedures perform operations on complex numbers.
#  
#       complex(r,i)    create complex number with real part r and
#                       imaginary part i
#  
#       cpxabs(z)       compute absolute value of complex number z
#
#       cpxadd(z1, z2)  add complex numbers z1 and z2
#
#       cpxconj(z)      compute conjugate of complex number z
#  
#       cpxdiv(z1, z2)  divide complex number z1 by complex number z2
#  
#       cpxmul(z1, z2)  multiply complex number z1 by complex number z2
#  
#       cpxsub(z1, z2)  subtract complex number z2 from complex number z1
#  
#       cpxstr(z)      convert complex number z to string representation
#  
#       strcpx(s)      convert string representation s of complex
#                      number to complex number
#  
############################################################################

record complex(rpart, ipart)

procedure strcpx(s)			#: convert string to complex number

   s ? {
      ="(" | fail
      return complex(numeric(upto('+-')),
         2(move(1), numeric(upto(')')), tab(-1)))
      }

end

procedure cpxstr(z)			#: return complex number as string

   if z.ipart < 0 then return "(" ||  z.rpart || z.ipart || "i)"
   else return "(" ||  z.rpart || "+" || z.ipart || "i)"

end

procedure cpxadd(z1, z2)		#: complex add

   return complex(z1.rpart + z2.rpart, z1.ipart + z2.ipart)

end

procedure cpxsub(z1, z2)		#: complex subtract

   return complex(z1.rpart - z2.rpart, z1.ipart - z2.ipart)

end

procedure cpxmul(z1, z2)		#: complex multiply

   return complex(z1.rpart * z2.rpart - z1.ipart * z2.ipart,
      z1.rpart * z2.ipart + z1.ipart * z2.rpart)

end

procedure cpxdiv(z1, z2)		#: complex divide
   local denom

   denom := z2.rpart ^ 2 + z2.ipart ^ 2

   return complex((z1.rpart * z2.rpart + z1.ipart * z2.ipart) / denom,
      (z1.ipart * z2.rpart - z1.rpart * z2.ipart) / denom)

end

procedure cpxconj(z)			#: complex conjugate

   return complex(z.rpart, -z.ipart)

end

procedure cpxabs(z)			#: complex absolute value

   return sqrt(z.rpart ^ 2 + z.ipart ^ 2)

end