131 lines
3.9 KiB
Plaintext
131 lines
3.9 KiB
Plaintext
'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)
|