'complex type and operators for it
type complex
    real as double
    imag as double
end type

operator + ( a as complex, b as complex ) as complex
    dim as complex r
    r.real = a.real + b.real
    r.imag = a.imag + b.imag
    return r
end operator

operator * ( a as complex, b as complex ) as complex
    dim as complex r
    r.real = a.real*b.real - a.imag*b.imag
    r.imag = a.real*b.imag + b.real*a.imag
    return r
end operator

operator = ( a as complex, b as complex ) as boolean
    if not a.real = b.real then return false
    if not a.imag = b.imag then return false
    return true
end operator

function complex_conjugate( a as complex ) as complex
    dim as complex r
    r.real =  a.real
    r.imag = -a.imag
    return r
end function

'matrix type and operations for it
'reuses code from the matrix multiplication task
type Matrix
    dim as complex m( any , any )
    declare constructor ( )
    declare constructor ( byval x as uinteger )
end type

constructor Matrix ( )
end constructor

constructor Matrix ( byval x as uinteger )
    redim this.m( x - 1 , x - 1 )
end constructor

operator * ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    redim ret.m( ubound( a.m , 1 ) , ubound( a.m , 1 ) )
    for i = 0 to ubound( a.m , 1 )
        for j = 0 to ubound( b.m , 2 )
            for k = 0 to ubound( b.m , 1 )
                ret.m( i , j ) += a.m( i , k ) * b.m( k , j )
            next k
        next j
    next i
    return ret
end operator

function conjugate_transpose( byref a as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j
    redim ret.m( ubound( a.m , 1 ) , ubound( a.m , 1 ) )
    for i = 0 to ubound( a.m , 1 )
        for j = 0 to ubound( a.m , 2 )
            ret.m( i, j ) = complex_conjugate(a.m( j, i ))
        next j
    next i
    return ret
end function

'tests if matrices are unitary, hermitian, or normal

operator = (byref a as Matrix, byref b as matrix) as boolean
    dim as integer i, j
    if ubound(a.m, 1) <> ubound(b.m, 1) then return false
    for i = 0 to ubound( a.m , 1 )
        for j = 0 to ubound( a.m , 2 )
            if not a.m(i,j)=b.m(i,j) then return false
        next j
    next i
    return true
end operator

function is_identity( byref a as Matrix ) as boolean
    dim as integer i, j
    for i = 0 to ubound( a.m , 1 )
        for j = 0 to ubound( a.m , 2 )
            if i  = j and ( not a.m(i,j).real = 1.0 or not a.m(i,j).imag = 0.0 ) then return false
            if i <> j and ( not a.m(i,j).real = 0.0 or not a.m(i,j).imag = 0.0 ) then return false
        next j
    next i
    return true
end function

function is_hermitian( byref a as Matrix ) as boolean
    if a = conjugate_transpose(a) then return true
    return false
end function

function is_normal( byref a as Matrix ) as boolean
    dim as Matrix aa = conjugate_transpose(a)
    if a*aa = aa*a then return true else return false
end function

function is_unitary( byref a as Matrix ) as boolean
    dim as Matrix aa = conjugate_transpose(a)
    if not is_identity( a*aa ) or not is_identity( aa*a ) then return false
    return true
end function

'''now some example matrices
dim as Matrix A = Matrix(2)   'an identity matrix
A.m(0,0).real = 1.0 : A.m(0,0).imag = 0.0 : A.m(0,1).real = 0.0 : A.m(0,1).imag = 0.0
A.m(1,0).real = 0.0 : A.m(1,0).imag = 0.0 : A.m(1,1).real = 1.0 : A.m(1,1).imag = 0.0

dim as Matrix B = Matrix(2)   'a hermitian matrix
B.m(0,0).real = 1.0 : B.m(0,0).imag = 0.0 : B.m(0,1).real = 1.0 : B.m(0,1).imag = -1.0
B.m(1,0).real = 1.0 : B.m(1,0).imag = 1.0 : B.m(1,1).real = 1.0 : B.m(1,1).imag = 0.0

dim as Matrix C = Matrix(2)   'a random matrix
C.m(0,0).real = rnd : C.m(0,0).imag = rnd : C.m(0,1).real = rnd : C.m(0,1).imag = rnd
C.m(1,0).real = rnd : C.m(1,0).imag = rnd : C.m(1,1).real = rnd : C.m(1,1).imag = rnd

print is_hermitian(A), is_normal(A), is_unitary(A)
print is_hermitian(B), is_normal(B), is_unitary(B)
print is_hermitian(C), is_normal(C), is_unitary(C)