FBComplex: complex number library
FBComplex: complex number library
http://www.unilim.fr/pages_perso/jean.d ... omplex.zip
This is a library for computing with complex numbers. It implements arithmetic operators (+, -, *, /, ^) and most transcendental functions (Log, Exp, Roots, Trigonometric, Hyperbolic, Gamma). It can be used with FBMath, or alone.
The archive includes:
- The source code of the library (fbcomplex.bas) and its header file (fbcomplex.bi)
- A precompiled library for Windows (libfbcomplex.a)
- A demo program (testcomp.bas)
- A PDF documentation
The library should be useful for creating fractal graphics, or for translating Fortran codes.
This is a library for computing with complex numbers. It implements arithmetic operators (+, -, *, /, ^) and most transcendental functions (Log, Exp, Roots, Trigonometric, Hyperbolic, Gamma). It can be used with FBMath, or alone.
The archive includes:
- The source code of the library (fbcomplex.bas) and its header file (fbcomplex.bi)
- A precompiled library for Windows (libfbcomplex.a)
- A demo program (testcomp.bas)
- A PDF documentation
The library should be useful for creating fractal graphics, or for translating Fortran codes.
Re: FBComplex: complex number library
Thanks jdebord.
This is easily setup (Windows) and works well.
10/10 or 10 stars.
I'm trying out the Gram-determinant method for Eigenvalues, If I need some functions I'll be able to pull them out of your source.
The gram-determinant method is analytical and needs no iteration.
I have been messing about with a dll from QB64, libquadmath-0.dll.
It's a gcc thing :
http://gcc.gnu.org/onlinedocs/libquadma ... -constants
I got the .def from tiny c, but as yet cannot get a runner.
By the way, when we awake in the morn (God willing), we may sing "here comes Summer", in the shower, for Winter stands behind us at 05.30 GMT.
Last year it was -15 centigrade, today was +10
Thanks again.
This is easily setup (Windows) and works well.
10/10 or 10 stars.
I'm trying out the Gram-determinant method for Eigenvalues, If I need some functions I'll be able to pull them out of your source.
The gram-determinant method is analytical and needs no iteration.
I have been messing about with a dll from QB64, libquadmath-0.dll.
It's a gcc thing :
http://gcc.gnu.org/onlinedocs/libquadma ... -constants
I got the .def from tiny c, but as yet cannot get a runner.
By the way, when we awake in the morn (God willing), we may sing "here comes Summer", in the shower, for Winter stands behind us at 05.30 GMT.
Last year it was -15 centigrade, today was +10
Thanks again.
Thank you Richard and dodicat for your comments !
@Richard:
I have now set all the implicit attributes, according to what I think is the default, i. e. :
- BYREF for UDT (Complex type) and STRING
- BYVAL for INTEGER and DOUBLE
@dodicat:
- Do you have the reference of this Gram-determinant method? I have never heard of it.
- I could not find libquadmat in my MinGW installation. I will check the QB64 distribution. If I understand correctly, there is no complex type in TCC.
- In France too the weather is much less cold than last year but it is rainy :)
@Richard:
I have now set all the implicit attributes, according to what I think is the default, i. e. :
- BYREF for UDT (Complex type) and STRING
- BYVAL for INTEGER and DOUBLE
@dodicat:
- Do you have the reference of this Gram-determinant method? I have never heard of it.
- I could not find libquadmat in my MinGW installation. I will check the QB64 distribution. If I understand correctly, there is no complex type in TCC.
- In France too the weather is much less cold than last year but it is rainy :)
For the Gram determinant ~ eigenvalues:
http://en.wikipedia.org/wiki/Gram_schmidt
(near the bottom).
I simply used tcc to extract the .def file to see the exports.
It's representing the type __float128 in fb which I find difficult.
The QB64 distribution has various other files with the dll.
Perhaps this dll is only a helper for 64 bit, I don't know.
Since the dll is GCC based, and fb is GCC inclined, I thought it would be a possibility to use this dll for Quad precision.
But we all know what thought did:
Planted a feather to grow a hen!
http://en.wikipedia.org/wiki/Gram_schmidt
(near the bottom).
I simply used tcc to extract the .def file to see the exports.
It's representing the type __float128 in fb which I find difficult.
The QB64 distribution has various other files with the dll.
Perhaps this dll is only a helper for 64 bit, I don't know.
Since the dll is GCC based, and fb is GCC inclined, I thought it would be a possibility to use this dll for Quad precision.
But we all know what thought did:
Planted a feather to grow a hen!
Sorry SARG, and I was born in 48, probably too young to pick up on it.SARG wrote:@Dodicat
Since 1948 we must use °Celius not °centigrade for temperature. It's a current error.
By the way, I like your last piece of code (Marcov/TJF), funny and very appropriate.
The last piece of code seems controversial.
@Richard
It was only a bit of fun, I'd had a couple of jars (Christmas and all that), nothing sinister.
@jdebord
May I be pedantic.
You've got complex^complex
-------------complex^ integer
-------------complex^double
----------but not real^complex
Just for the sake of Euler's identity, but if e is complex it's ok.
Code: Select all
#INCLUDE "fbcomplex.bi"
dim as complex e=(exp(1),0)
dim as complex i=(0,1)
dim as double pi=4*atn(1)
cprint (e^(i*pi)+1,"###.###")
sleep
I have added the real^complex operator but the web site is not updated yet.
I wanted to add functions for polynomials and rational fractions but I don't know how the arrays are passed. For instance the following function:
compiles fine, but when I add either BYVAL or BYREF before Coef() I get an error :
I thought that arrays are passed by reference.
I wanted to add functions for polynomials and rational fractions but I don't know how the arrays are passed. For instance the following function:
Code: Select all
FUNCTION CPoly(BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg AS INTEGER) AS Complex
' Evaluates the polynomial :
' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg
DIM AS INTEGER I
DIM AS Complex P
P = Coef(Deg)
FOR I = Deg - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
RETURN P
END FUNCTION
Code: Select all
Illegal specification, at parameter 2 (Coef) of CPoly()
Extract of manual about keyword 'Function':jdebord wrote:I thought that arrays are passed by reference.
...
Passing Arguments : ... Array parameters are specified by following an identifier with an empty parenthesis. Note that array parameters are always ByRef and the ByRef keyword is neither required nor allowed for array parameters. When calling a function with an array argument the parenthesis must be supplied there too...
Hi jdebordjdebord wrote:I have added the real^complex operator but the web site is not updated yet.
I wanted to add functions for polynomials and rational fractions but I don't know how the arrays are passed. For instance the following function:
compiles fine, but when I add either BYVAL or BYREF before Coef() I get an error :Code: Select all
FUNCTION CPoly(BYREF Z AS Complex, _ Coef() AS Complex, _ BYVAL Deg AS INTEGER) AS Complex ' Evaluates the polynomial : ' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg DIM AS INTEGER I DIM AS Complex P P = Coef(Deg) FOR I = Deg - 1 TO 0 STEP -1 P = P * Z + Coef(I) NEXT I RETURN P END FUNCTION
I thought that arrays are passed by reference.Code: Select all
Illegal specification, at parameter 2 (Coef) of CPoly()
The way you have it, -w pedantic gives no error.
So is this not just correct.
For my own Cpoly , to save a parameter I get deg as ubound(the input array).
Although yours is more compact and probaly faster.
-w pedantic gives no errors in this comparison
Code: Select all
type complex
as double re,im
end type
operator +(byval n1 As complex,byval n2 As complex) As complex
Return type<complex>(n1.re+n2.re,n1.im+n2.im)'n
End operator
operator *(byval n1 As complex,byval n2 As complex) As complex
Return type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)'n
End operator
Function Cpoly2(coff() As complex,byval number As complex)As complex
'evaluates the polynomial
Dim pol As complex
Dim deg As Integer=Ubound(coff)
For count as integer = 1 To DEG + 1
var temp=type<complex>(1,0)
for k as integer=1 to count-1
temp=temp*number
next k
pol = pol + coff(count-1) * temp
Next count
Cpoly2 = pol
End Function
FUNCTION CPoly(BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg AS INTEGER) AS Complex
' Evaluates the polynomial :
' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg
DIM AS INTEGER I
DIM AS Complex P
P = Coef(Deg)
FOR I = Deg - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
RETURN P
END FUNCTION
dim as complex poly(4)
dim as complex answer,pt
for z as integer=1 to 10
randomize
poly(0)=type<complex>(rnd*3-rnd*3,rnd*.08-rnd*.08)
poly(1)=type<complex>(-5,.2)
poly(2)=type<complex>(0,-8)
poly(3)=type<complex>(.5,0)
poly(4)=type<complex>(rnd*5-rnd*5,rnd*.5-rnd*.5)
pt=type<complex>(rnd*1-rnd*1,rnd*2-rnd*2)
answer=cpoly(pt,poly(),4)
print "jdb ";answer.re,answer.im
answer=cpoly2(poly(),pt)
print "dod ";answer.re,answer.im
print
next z
sleep
Hi jdebord
A little Christmas bundle.
I got some help from another Frenchman
http://jean-pierre.moreau.pagesperso-or ... basic.html
Set up for polynomial roots (real or complex polynomial).
And eigenvalues/vectors.
I have used your Cpoly, and my operators, but easily changed to your operators, 15 minutes maybe.
All the main stuff is about 400 lines, the rest of the lines are just for the examples.
To input a complex component just key in e.g. 9,-5
It's maybe a bit untidy, do as you wish with it, after all, the essence came from France, and I am just sending it back across the Channel.
Merry Christmas.
A little Christmas bundle.
I got some help from another Frenchman
http://jean-pierre.moreau.pagesperso-or ... basic.html
Set up for polynomial roots (real or complex polynomial).
And eigenvalues/vectors.
I have used your Cpoly, and my operators, but easily changed to your operators, 15 minutes maybe.
All the main stuff is about 400 lines, the rest of the lines are just for the examples.
To input a complex component just key in e.g. 9,-5
It's maybe a bit untidy, do as you wish with it, after all, the essence came from France, and I am just sending it back across the Channel.
Merry Christmas.
Code: Select all
Type complex
As Double re,im
End Type
Operator +(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re+n2.re,n1.im+n2.im)
End Operator
Operator -(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re-n2.re,n1.im-n2.im)
End Operator
Operator *(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)
End Operator
Operator /(n1 As complex,n2 As complex) As complex
Dim As Double d = n2.re*n2.re+n2.im*n2.im
Return Type<complex>((n1.re*n2.re+n1.im*n2.im)/d,(n1.im*n2.re - n1.re*n2.im)/d)
End Operator
Declare Sub geteigens(m() As complex,evals() As complex,evects() As complex,Byref IER As Integer)
Declare Sub COMPLEXEIG (A() As Double,Z() As Double,T() As Double,U() As Double,N As Integer, NN As Integer, IER As Integer)
'_____________ All above this line is necessary to get eigenvalues/vectors ________
'_______________ These for polyroots __________________________
declare sub prepare_polynomial(poly()as complex)
Declare Sub makecompanion(polynomium() As complex , mat() As complex)
'These Four subs below are for examples and printouts only
Declare Function Cpoly(coff() As complex,number As complex)As complex 'evaluate poly(YOURS)
Declare Sub cprint(x As complex,places As Integer=10) 'print a rounded complex number
Declare Function decround(a As Double , b As Integer) As Double 'the rounder function
declare Sub matmultvect(m1() As complex,m2() As complex,ans() As complex) 'MATRIX x VECTOR
'__________________ EXAMPLE ROOTS ____________________________
Dim deg As Integer
do
Input "Enter polynomial degree (for a random polynomial) ",deg
loop until deg>0
Print
'print "Input coefficients thus: 6+7i is entered as 6,7"
'print "If coefficient is real then just enter as is, i.e. .55 is just entered as .55"
print
Dim As complex poly(deg),copypoly(deg)
start:
For z As Integer=deg To 0 Step -1
If z>0 Then
'' Print "Enter coefficient of x^";Str(z);
''Input " ",poly(z).re,poly(z).im
poly(z)=type<complex>(rnd*2,rnd*3)
copypoly(z)=poly(z)
if poly(deg).re=0 and poly(deg).im = 0 then
''print "First coefficient must not be zero"
goto start
end if
End If
If z=0 Then
''Print "Enter constant term ";
''Input " ",poly(z).re,poly(z).im
poly(z)=type<complex>(rnd*2,rnd*3)
copypoly(z)=poly(z)
End If
Next z
'THE WORKING VARIABLES
Dim As complex companion(0,0)
Dim As Integer IER 'ERROR INDICATOR
Redim As complex values(0),vectors(0,0) 'Eigenvalues/vectors
prepare_polynomial(poly())
makecompanion poly(),companion()
'The eigenvalues of the companion matrix are the roots
geteigens(companion(),values(),vectors(),IER) 'GET EIGENVALUES OF COMPANION
'(DONT USE VECTORS FOR POLY ROOTS)
Print
Print "ROOTS (To 15 places):"
print
For I As Integer=1 To Ubound (values)
cprint(values(I),15)
Next i
Print
Print " Error code (must be zero): "; IER
Print
Print "POLYNOMIAL VALUES AT ROOTS (to six places), should be (0, 0i) for each one"
print
For z As Integer=1 To deg
cprint(cpoly(copypoly(),values(z)),6)
Next z
print
'EXAMPLE EIGENVALUES/EIGENVECTORS
print "Eigenvectors example (KEY inputs, so not too big!)"
print
reDim As complex matrix(0,0),eigenvalues(0),eigenvectors(0,0)
Redim As complex product(0) 'for checking products
input "Enter dimension for key input",deg
reDim As complex matrix(1 to deg,1 to deg)
redim eigenvalues(1 to deg)
redim eigenvectors(1 to deg,1 to deg)
redim product(1 to deg)
For _row As Integer=1 To deg
For _column As Integer=1 To deg
Print "Enter (";_row;",";_column;")";
Input matrix(_row,_column).re,matrix(_row,_column).im
'matrix(_row,_column).re=Rnd*7
'matrix(_row,_column).im=Rnd*8
Next _column
Print
Next _row
Print "MATRIX IS:"
For aa As Integer=1 To deg
For bb As Integer=1 To deg
Print "(";Str(decround(matrix(aa,bb).re,10));",";Str(decround(matrix(aa,bb).im,10));"i) ";
Next bb
Print
Next aa
Print
geteigens(matrix(),eigenvalues(),eigenvectors(),IER)
Print "Eigenvalues:"
For z As Integer=1 To deg
cprint eigenvalues(z)
Next z
Print
Print
'display eigenvectors
Print" Corresponding Eigenvectors:"
For z2 As Integer=1 To deg
For z As Integer=1 To deg
cprint eigenvectors(z,z2)
Next z
Print
Next z2
'FOR SCREEN DISPLAY ONLY
print "CHECK"
print
Print "Matrix x eigenvector Eigenvalue x eigenvector"
print
Dim As complex temp(1 To deg)
Dim As complex temp2
redim as string MxE(0),ExE(0) 'to hold values for screen display
dim as integer count
'multiply eigenvector by matrix and save in array MxE() of string
For z2 As Integer=1 To deg
For z As Integer=1 To deg
temp(z)=eigenvectors(z,z2)
Next z
matmultvect(matrix(),temp(),product())'get the product matrix (array)
For z3 As Integer=1 To deg
count=count+1
redim preserve MxE(count) 'save array for display
MxE(count)="("+str(decround(product(z3).re,8))+","+str(decround(product(z3).im,8))+"i)"
Next z3
count=count+1
redim preserve MxE(count)
MxE(count)=" "
Next z2
'multiply each eigenvector element by its corresponding eigenvalue and save in ExE()
count=0
For z2 As Integer=1 To deg
For z As Integer=1 To deg
temp2=eigenvalues(z2)*eigenvectors(z,z2)'mult each eigenvector component
count=count+1
redim preserve ExE(count)'save here for display
ExE(count)="("+str(decround(temp2.re,8))+","+str(decround(temp2.im,8))+"i)"
Next z
count=count+1
redim preserve ExE(count)
ExE(count)=" "
Next z2
'print comparison on screen
for z as integer=1 to count
print MxE(z);tab(40);ExE(z)
next z
Print "Error = ";IER
print "DONE"
Sleep
End 'of examples
'______________
'____________________________________________
sub prepare_polynomial(poly()as complex)
dim deg as integer=ubound(poly)
Dim As complex lead=poly(deg) 'LEADING COEFFICIENT
For z As Integer=deg To 0 Step -1
poly(z)=poly(z)/lead 'PREPARE POLY FOR COMPANION MATRIX
Next z
end sub
Sub geteigens(m() As complex,evals() As complex,evects() As complex,Byref IER As Integer)
Dim As Integer n=Ubound(m),I,J
Dim As Double A(n,n),Z(n,n),T(n,n),U(n,n)
For z1 As Integer=1 To n
For z2 As Integer=1 To n
a(z1,z2)=m(z1,z2).re
z(z1,z2)=m(z1,z2).im
Next z2
Next z1
ComplexEIG A(),Z(),T(),U(),N, 25, IER
Redim evals(n)
For I = 1 To N
evals(I).re=A(I,I)
evals(I).im=Z(I,I)
Next I
Redim evects(n,n)
For I=1 To n
For J=1 To n
evects(i,j).re=T(i,j)
evects(i,j).im=U(i,j)
Next j
Next i
End Sub
Sub COMPlexEIG (A() As Double,Z() As Double,T() As Double,U() As Double,N As Integer, NN As Integer, IER As Integer)
#macro complexabs(r,i)
Sqr(r*r+i*i)
#endmacro
#macro complexdiv2(ar,ai,br,bi,zr,zi)
Scope
var YR = BR
var YI = BI
If Abs(YR) > Abs(YI) Then
var W = YI / YR
YR = W * YI + YR
ZR = (AR + W * AI) / YR
ZI = (AI - W * AR) / YR
Else
var W = YR / YI
YI = W * YR + YI
ZR = (W * AR + AI) / YI
ZI = (W * AI - AR) / YI
End If
End Scope
#endmacro
Dim As Double EPS,EPS1,TAU,W1,EMAX,W2,G,HR,HI,HJ,T1,T2,BI,ER,EI,DR,DI
Dim As Double SW,C,S,D,DE,ROOT2,ROOT1,SIG,CA,SA,SX,CX,DN,B,E
Dim As Double COT2X,COTX,ETA,TSE,COS2A,SIN2A,ROOT,CB,CH,SB,SH,CN,TANH
Dim As Double C1R,C2R,C1I,C2I,S1R,S1I,S2R,S2I,AKI,AMI,ZKI,ZMI,AIK,AIM,ZIK,ZIM,TIK,TIM,UIK,UIM
Dim As Double VR,VI,ZM,ZI,IM,TR,TI,BR
Dim As Integer MARK,I,J,K,IT,M
Dim As Double EN(2*N)
EPS = 1#
ten: EPS = .5 * EPS
EPS1 = EPS + 1#
If EPS1 > 1# Then Goto ten
MARK = 0 'FALSE
' INITIALIZE EIGENVECTORS
For I = 1 To N
T(I, I) = 1#
U(I, I) = 0#
For J = I + 1 To N
T(I, J) = 0#
T(J, I) = 0#
U(I, J) = 0#
U(J, I) = 0#
Next J
Next I
IT = 0
Twenty: IT = IT + 1
' SAFETY TEST IN CASE OF NO CONVERGENCE
If IT > NN Then Goto Ninety
If MARK = 1 Then Goto NinetyFive '1 = TRUE
' DEFINE CONVERGENCE CRITERIUM
TAU = 0#
For K = 1 To N
W1 = 0#
For I = 1 To N
If I <> K Then W1 = W1 + Abs(A(I, K)) + Abs(Z(I, K))
Next I
TAU = TAU + W1
EN(K) = W1 + Abs(A(K, K)) + Abs(Z(K, K))
Next K
' PERMUTE LINES AND COLUMNS
For K = 1 To N - 1
EMAX = EN(K)
I = K
For J = K + 1 To N
If EN(J) > EMAX Then
EMAX = EN(J)
I = J
End If
Next J
If I <> K Then
EN(I) = EN(K)
For J = 1 To N
W2 = A(K, J)
A(K, J) = A(I, J)
A(I, J) = W2
W2 = Z(K, J)
Z(K, J) = Z(I, J)
Z(I, J) = W2
Next J
For J = 1 To N
W2 = A(J, K)
A(J, K) = A(J, I)
A(J, I) = W2
W2 = Z(J, K)
Z(J, K) = Z(J, I)
Z(J, I) = W2
W2 = T(J, K)
T(J, K) = T(J, I)
T(J, I) = W2
W2 = U(J, K)
U(J, K) = U(J, I)
U(J, I) = W2
Next J
End If
Next K
' CONVERGENCE IF TAU < 100*EPS
If TAU < 100# * EPS Then Goto NinetyFive '100
' BEGIN ITERATIONS
MARK = 1 'TRUE
For K = 1 To N - 1
For M = K + 1 To N
G = 0#
HR = 0#
HJ = 0#
HI = 0#
For I = 1 To N
If I <> K And I <> M Then
HR = HR + A(K, I) * A(M, I) + Z(K, I) * Z(M, I) - A(I, K) * A(I, M) - Z(I, K) * Z(I, M)
HI = HI + Z(K, I) * A(M, I) - A(K, I) * Z(M, I) - A(I, K) * Z(I, M) + Z(I, K) * A(I, M)
T1 = A(I, K) * A(I, K) + Z(I, K) * Z(I, K) + A(M, I) * A(M, I) + Z(M, I) * Z(M, I)
T2 = A(I, M) * A(I, M) + Z(I, M) * Z(I, M) + A(K, I) * A(K, I) + Z(K, I) * Z(K, I)
G = G + T1 + T2
HJ = HJ - T1 + T2
End If
Next I
BR = A(K, M) + A(M, K)
BI = Z(K, M) + Z(M, K)
ER = A(K, M) - A(M, K)
EI = Z(K, M) - Z(M, K)
DR = A(K, K) - A(M, M)
DI = Z(K, K) - Z(M, M)
T1 = BR * BR + EI * EI + DR * DR
T2 = BI * BI + ER * ER + DI * DI
If T1 >= T2 Then
SW = 1#
C = BR
S = EI
D = DR
DE = DI
ROOT2 = Sqr(T1)
Else
SW = -1#
C = BI
S = -ER
D = DI
DE = DR
ROOT2 = Sqr(T2)
End If
ROOT1 = Sqr(S * S + C * C)
SIG = 1#
If D < 0# Then SIG = -1#
CA = 1#
If C < 0# Then CA = -1#
SA = 0#
If ROOT1 < EPS Then
SX = 0#
SA = 0#
CX = 1#
CA = 1#
If SW > 0# Then
E = ER
B = BI
Else
E = EI
B = -BR
End If
DN = D * D + DE * DE
Goto Sixtyfive
End If
If Abs(S) > EPS Then
CA = C / ROOT1
SA = S / ROOT1
End If
COT2X = D / ROOT1
COTX = COT2X + SIG * Sqr(1# + COT2X * COT2X)
SX = SIG / Sqr(1# + COTX * COTX)
CX = SX * COTX
ETA = (ER * BR + BI * EI) / ROOT1
TSE = (BR * BI - ER * EI) / ROOT1
T1 = SIG * (TSE * D - ROOT1 * DE) / ROOT2
T2 = (D * DE + ROOT1 * TSE) / ROOT2
DN = ROOT2 * ROOT2 + T2 * T2
T2 = HJ * CX * SX
COS2A = CA * CA - SA * SA
SIN2A = 2! * CA * SA
W1 = HR * COS2A + HI * SIN2A
W2 = HI * COS2A - HR * SIN2A
HR = CX * CX * HR - SX * SX * W1 - CA * T2
HI = CX * CX * HI + SX * SX * W2 - SA * T2
B = SW * T1 * CA + ETA * SA
E = CA * ETA - SW * T1 * SA
' ROOT1 < EPS
Sixtyfive: S = HR - SIG * ROOT2 * E
C = HI - SIG * ROOT2 * B
ROOT = Sqr(C * C + S * S)
If ROOT < EPS Then
CB = 1#
CH = 1#
SB = 0#
SH = 0#
Goto Seventy
End If
CB = -C / ROOT
SB = S / ROOT
T2 = CB * B - E * SB
CN = T2 * T2
TANH = ROOT / (G + 2# * (CN + DN))
CH = 1# / Sqr(1# - TANH * TANH)
SH = CH * TANH
' ROOT < EPS
Seventy: W1 = SX * SH * (SA * CB - SB * CA)
C1R = CX * CH - W1
C2R = CX * CH + W1
C1I = -SX * SH * (CA * CB + SA * SB)
C2I = C1I
W2 = SX * CH * CA
W1 = CX * SH * SB
S1R = W2 - W1
S2R = -W2 - W1
W2 = SX * CH * SA
W1 = CX * SH * CB
S1I = W2 + W1
S2I = W2 - W1
W1 = Sqr(S1R * S1R + S1I * S1I)
W2 = Sqr(S2R * S2R + S2I * S2I)
If W1 > EPS Or W2 > EPS Then
' BEGIN TRANSFORMATIONS
MARK = 0 'FALSE
For I = 1 To N
AKI = A(K, I)
AMI = A(M, I)
ZKI = Z(K, I)
ZMI = Z(M, I)
A(K, I) = C1R * AKI - C1I * ZKI + S1R * AMI - S1I * ZMI
Z(K, I) = C1R * ZKI + C1I * AKI + S1R * ZMI + S1I * AMI
A(M, I) = S2R * AKI - S2I * ZKI + C2R * AMI - C2I * ZMI
Z(M, I) = S2R * ZKI + S2I * AKI + C2R * ZMI + C2I * AMI
Next I
For I = 1 To N
AIK = A(I, K)
AIM = A(I, M)
ZIK = Z(I, K)
ZIM = Z(I, M)
TIK = T(I, K)
TIM = T(I, M)
UIK = U(I, K)
UIM = U(I, M)
A(I, K) = C2R * AIK - C2I * ZIK - S2R * AIM + S2I * ZIM
Z(I, K) = C2R * ZIK + C2I * AIK - S2R * ZIM - S2I * AIM
A(I, M) = C1R * AIM - C1I * ZIM - S1R * AIK + S1I * ZIK
Z(I, M) = C1R * ZIM + C1I * AIM - S1R * ZIK - S1I * AIK
T(I, K) = C2R * TIK - C2I * UIK - S2R * TIM + S2I * UIM
U(I, K) = C2R * UIK + C2I * TIK - S2R * UIM - S2I * TIM
T(I, M) = C1R * TIM - C1I * UIM - S1R * TIK + S1I * UIK
U(I, M) = C1R * UIM + C1I * TIM - S1R * UIK - S1I * TIK
Next I
End If
' END TRANSFORMATIONS
Next M
Next K
' GO TO NEXT ITERATION
Goto Twenty
' NO CONVERGENCE '
Ninety: IER = 1
Exit Sub
' CONVERGENCE OK
NinetyFive: IER = 0
' SORT SOLUTIONS IN INCREASING ORDER
For J = 2 To N
VR = A(J, J)
VI = Z(J, J)
For K = 1 To N
EN(K) = T(K, J)
EN(K + N) = U(K, J)
Next K
For I = J - 1 To 1 Step -1
If ComplexABS(A(I, I), Z(I, I)) <= ComplexABS(VR, VI) Then Goto Ninetyseven
A(I + 1, I + 1) = A(I, I)
Z(I + 1, I + 1) = Z(I, I)
For K = 1 To N
T(K, I + 1) = T(K, I)
U(K, I + 1) = U(K, I)
Next K
Next I
I = 0
Ninetyseven: A(I + 1, I + 1) = VR
Z(I + 1, I + 1) = VI
For K = 1 To N
T(K, I + 1) = EN(K)
U(K, I + 1) = EN(K + N)
Next K
Next J
' NORMALIZE VECTORS (BIGGEST COMPONENT TO UNITY)
'exit sub' if you don't want to normalize
For J = 1 To N
ZM = 0#
For I = 1 To N
ZI = Abs(T(I, J)) + Abs(U(I, J))
If ZI >= ZM Then
IM = I
ZM = ZI
End If
Next I
ZM = T(IM, J)
ZI = U(IM, J)
For I = 1 To N
ComplexDIV2 (T(I, J), U(I, J), ZM, ZI, TR, TI)
T(I, J) = TR
U(I, J) = TI
Next I
Next J
End Sub
Sub makecompanion(polynomium() As complex , mat() As complex)
Dim As Integer di = Ubound(polynomium)'+1
Dim As complex one , _null
one.re = 1
Redim mat(1 To di , 1 To di) As complex
For a As Integer = 1 To di
mat(1 , di) = _null - polynomium(0)
For b As Integer = 1 To di
If a = b + 1 Then mat(a , b) = one
If b > 1 Then
mat(b , di) = _null - polynomium(b - 1)
End If
Next b
Next a
End Sub
'__________________________________________________________________
'THE HELPER SUBS FOR EXAMPLES
FUNCTION CPoly(coef() as complex,BYREF Z AS Complex) as complex
dim as integer deg=ubound(coef)
' Evaluates the polynomial :
' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg
DIM AS Complex P
P = Coef(Deg)
FOR I as integer= Deg - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
RETURN P
END FUNCTION
Function decround(a As Double , b As Integer) As Double
Dim y As Double
Dim i As Double
Dim r As Double
y = (Abs(a) - Int(Abs(a))) * (10^b)
i = Int(y)
y = y - i
If y >= .5 Then i = i + 1
i = i / (10^b)
r = Int(Abs(a)) + i
If a < 0 Then r = - r
decround = r
End Function
Sub cprint(x As complex,places as integer=10)
Print "(";decround(x.re,places),decround(x.im,places);"i)"
'Print "(";x.re,x.im;"i)"
End Sub
Sub matmultvect(m1() As complex,m2() As complex,ans() As complex) 'MATRIX x VECTOR
Dim s As complex
Dim n As Double=Ubound(m1)
Redim ans(1 To n)
For i As Integer=1 To n
s.re=0:s.im=0
For k As Integer = 1 To n
s=s+m1(i,k)*m2(k)
Next k
ans(i)=s
Next i
End Sub
Thanks dodicat, and merry Christmas !
In FBMath there are procedures for solving equations up to degree 4 with analytic formulas, and use the companion matrix for higher degrees.
J. P. Moreau's code for eigenvalues looks similar to the EISPACK code which I used for the eigenvalues of real matrices.
I still have to adapt the procedures to complex matrices.
In the meantime, here is the rational fraction:
In FBMath there are procedures for solving equations up to degree 4 with analytic formulas, and use the companion matrix for higher degrees.
J. P. Moreau's code for eigenvalues looks similar to the EISPACK code which I used for the eigenvalues of real matrices.
I still have to adapt the procedures to complex matrices.
In the meantime, here is the rational fraction:
Code: Select all
FUNCTION CFrac(BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg1 AS INTEGER, _
BYVAL Deg2 AS INTEGER) AS Complex
' Evaluates the rational fraction :
'
' Coef(0) + Coef(1) * Z + ... + Coef(Deg1) * Z^Deg1
' F(Z) = -----------------------------------------------------
' 1 + Coef(Deg1+1) * Z + ... + Coef(Deg1+Deg2) * Z^Deg2
DIM AS INTEGER I
DIM AS Complex P = Coef(Deg1), Q = (0, 0)
FOR I = Deg1 - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
Q = (Q + Coef(I)) * Z
NEXT I
RETURN P / (1 + Q)
END FUNCTION
Here is an improved version, with both real and complex coefficients:
Code: Select all
DECLARE FUNCTION CPoly OVERLOAD (BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg AS INTEGER) AS Complex
DECLARE FUNCTION CPoly (BYREF Z AS Complex, _
Coef() AS DOUBLE, _
BYVAL Deg AS INTEGER) AS Complex
DECLARE FUNCTION CFrac OVERLOAD (BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg1 AS INTEGER, _
BYVAL Deg2 AS INTEGER) AS Complex
DECLARE FUNCTION CFrac (BYREF Z AS Complex, _
Coef() AS DOUBLE, _
BYVAL Deg1 AS INTEGER, _
BYVAL Deg2 AS INTEGER) AS Complex
' ------------------------------------------------------------------
FUNCTION CPoly(BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg AS INTEGER) AS Complex
' Polynomial with complex coefficients
' P(Z) = Coef(0) + Coef(1) * Z + Coef(2) * Z^2 +... + Coef(Deg) * Z^Deg
DIM AS INTEGER I
DIM AS Complex P
P = Coef(Deg)
FOR I = Deg - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
RETURN P
END FUNCTION
FUNCTION CPoly(BYREF Z AS Complex, _
Coef() AS DOUBLE, _
BYVAL Deg AS INTEGER) AS Complex
' Polynomial with real coefficients
DIM AS INTEGER I
DIM AS Complex P = Cmplx(Coef(Deg), 0)
FOR I = Deg - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
RETURN P
END FUNCTION
FUNCTION CFrac(BYREF Z AS Complex, _
Coef() AS Complex, _
BYVAL Deg1 AS INTEGER, _
BYVAL Deg2 AS INTEGER) AS Complex
' Rational fraction with complex coefficients
'
' Coef(0) + Coef(1) * Z + ... + Coef(Deg1) * Z^Deg1
' F(Z) = --------------------------------------------------------------------
' Coef(Deg1 + 1) + Coef(Deg1+2) * Z + ... + Coef(Deg1+Deg2+1) * Z^Deg2
DIM AS INTEGER I
DIM AS Complex P = Coef(Deg1)
DIM AS Complex Q = Coef(Deg1 + Deg2 + 1)
FOR I = Deg1 - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
Q = Q * Z + Coef(I)
NEXT I
RETURN P / Q
END FUNCTION
FUNCTION CFrac(BYREF Z AS Complex, _
Coef() AS DOUBLE, _
BYVAL Deg1 AS INTEGER, _
BYVAL Deg2 AS INTEGER) AS Complex
' Rational fraction with real coefficients
DIM AS INTEGER I
DIM AS Complex P = Cmplx(Coef(Deg1), 0)
DIM AS Complex Q = Cmplx(Coef(Deg1 + Deg2 + 1), 0)
FOR I = Deg1 - 1 TO 0 STEP -1
P = P * Z + Coef(I)
NEXT I
FOR I = Deg1 + Deg2 TO Deg1 + 1 STEP -1
Q = Q * Z + Coef(I)
NEXT I
RETURN P / Q
END FUNCTION
Thanks jdebord.
Here's my complex matrix inverse/determinant taylored to your latest library.
I notice that to return the product of determinant and determinant of inverse gives 1.000 - 0.000*i by your Cprint, which is correct of course, but why the minus sign?
I've tested with -exx -w pedantic with no comments.
100 x 100 in a heartbeat.
If it's of any use then do as you wish with it.
Full pivoting on reals, so should be accurate.
I've multiplied matrix x inverse as a second check.
Here's my complex matrix inverse/determinant taylored to your latest library.
I notice that to return the product of determinant and determinant of inverse gives 1.000 - 0.000*i by your Cprint, which is correct of course, but why the minus sign?
I've tested with -exx -w pedantic with no comments.
100 x 100 in a heartbeat.
If it's of any use then do as you wish with it.
Full pivoting on reals, so should be accurate.
I've multiplied matrix x inverse as a second check.
Code: Select all
#INCLUDE "fbcomplex.bi"
'COMPLEX MATRIX INVERSE/DETERMINANT
declare Function DETERMINANT_INVERT(matrix() As complex,inverse() As complex) As complex ' determinant+inverse
declare Sub matrixmult(m1() As complex,m2() As complex,ans() As complex) ' multiply
declare Function decround (byval a As Double,byval b As Integer) As Double 'a decimal rounder
declare Function rr(byval first As Double, byval last As Double) As Double 'random generator
' ***************** EXAMPLE ************************
dim m as integer
start:
print "MATRIX WILL BE RANDOMLY FILLED"
Input "ENTER DIMENSION OF MATRIX or zero to quit) ";m
If m=0 Then End
Print
print "Please wait, calculating determinant, inverse and determinant of inverse"
Print
Redim matrix(1 To m,1 To m) As complex 'main matrix
Redim inverse(0,0) As complex 'THE FUNCTION DETERMINANT_INVERT() WILL FILL THIS
'MATRIX INPUT (randoms or manual)_________________________________
For aa As Integer=1 To m
For bb As Integer=1 To m
' Print "Enter (";aa;",";bb;")";
'Input matrix(aa,bb)
matrix(aa,bb).x=rr(-1,1)
matrix(aa,bb).y=rr(-1,1)
Next bb
'Print
Next aa
'_____________________________________________
Print
'DETERMINANT OUTPUT
Dim As complex d1,d2
d1=DETERMINANT_INVERT (matrix(),inverse()) 'returns determinant and inverse
Print "Determinant of matrix = ";d1.x;" ";d1.y;" i"
redim as complex dummy(0,0) 'only needed for array parameter
'Get determinant of the inverse
d2=DETERMINANT_INVERT(inverse(),dummy()) 'use the function again
Print "Determinant of inverse = ";d2.x;" ";d2.y;" i"
print
dim as complex d1Xd2
d1Xd2=d1*d2 'MULTIPLY DETERMINANTS
print "To 12 decimal places:"
print "Product of determinants:"
print
print "Product of determinants ";decround(d1Xd2.x,12);" ";decround(d1Xd2.y,12);" i"
'cprint(d1Xd2,"#####.############")
print
print "Now multiply matrix by it's inverse - press a key"
sleep
print
'INVERSE OUTPUT if required________________________________
'Print "INVERSE"
'For aa As Integer=1 To m
'For bb As Integer=1 To m
'Print "(";aa;",";bb;") = ";inverse(aa,bb).x;" "; inverse(aa,bb).y
'Next bb
'Print
'Next aa
'___________________________________________________
'NOW MULTIPLY MATRIX BY IT'S INVERSE TO CHECK
Redim product(0,0) As complex
matrixmult matrix(),inverse(),product()
Print "CHECK, MATRIX X INVERSE - SHOULD BE A UNIT MATRIX"
Print "Results for real part to 12 decimal places"
Print
For aa As Integer=1 To m
For bb As Integer=1 To m
Print decround(product(aa,bb).x,12); 'round to 12 places
Next bb
Print
Next aa
Print
Goto start
' ********************** END EXAMPLES *************************
'returns determinant and inverse of a complex matrix()
Function DETERMINANT_INVERT(matrix() As complex,inverse() As complex) As complex
'_____________________________________________________________
#macro setup()
dim as complex one=type<complex>(1,0),zero,mm,sign=type<complex>(1,0)
Dim As Long x2,y2 'counters for invert
Dim As Long x,y,z 'counters for lu
Dim det As complex=one 'determinant starter
Dim m As Long=Ubound(matrix)
Dim As Integer x1,k1,g,dimcount 'some counters
Redim h(1 To m) As complex 'column of a unit matrix
Redim l(1 To m,1 To m) As complex 'lower triangular
Redim b(1 To m,1 To m) As complex 'upper triangular
Redim iv(1 To m) As complex 'intermediate vector
Redim inverse(1 To m,1 To m) 'the inverse
Redim h2(1 To m,1 To m) As complex 'unit matrix
Redim solution(1 To m) As complex 'a solution column
For i As Long=1 To m
For j As Long=1 To m
b(i,j)=matrix(i,j) 'take a copy of the matrix
If i=j Then
l(i,j)=one 'make two unit matrices
h2(i,j)=one
Else
l(i,j)=zero
h2(i,j)=zero
Endif
Next j
Next i
#endmacro
'Everything setup
'____________________________________________________________
#macro cleanup()
Erase l,b,h,h2,iv,solution
#endmacro
' Annihilate arrays
'_____________________________________________________________
#macro lu()
For y=1 To m-1
If (b(y,y).x =0 and b(y,y).y=0) Then
pivot()
If (b(y,y).x and b(y,y).y=0) Then '
cleanup()
Return zero
End If
Endif
For x=y To m-1
l(x+1,y)=b(x+1,y)/b(y,y) 'l() is lower triangular matrix
mm=l(x+1,y)
b(x+1,y)=zero
For z=y+1 To m
b(x+1,z)=b(x+1,z)-mm*b(y,z) 'b() is upper triangular matrix
Next z
Next x
Next y
#endmacro
'Upper and Lower triangular matrices made up
'___________________________________________________________
#macro pivot()
For dimcount=0 To Ubound(matrix)-1
For x1=1 To Ubound(matrix)
For k1=1 To Ubound(matrix)-1
If k1+1+dimcount >ubound(matrix) Then Exit For
If x1+dimcount>UBound(matrix) Then Exit For
If Abs(b(k1+1+dimcount,1+dimcount).x)>Abs(b(x1+dimcount,1+dimcount).x) Then
'sign changes with each swap
sign=zero-sign 'equivalent to sign=-1*sign
For g As Long=1 To Ubound(matrix)
Swap b(k1+1+dimcount,g),b(x1+dimcount,g) 'swap lower triangle row
Swap h2(g,k1+1+dimcount),h2(g,x1+dimcount) 'swap unit row
Next g
End If
Next k1
Next x1
Next dimcount
#endmacro
'Biggest to top, smallest to the bottom, single rank sized reals
'__________________________________________________________
#macro ivector()
For n As Long=1 To m
iv(n)=h(n)
For j As Long=1 To n-1
iv(n)=iv(n)-l(n,j)*iv(j)
Next j
Next n
#endmacro
'intermediate vector made up
'___________________________________________________________
#macro fvector()
For n As Long=m To 1 Step -1
solution(n)=iv(n)/b(n,n)
For j As Long = m To n+1 Step -1
solution(n)=solution(n)-(b(n,j)*solution(j)/b(n,n))
Next j
Next n
#endmacro
'solution vector made up
'___________________________________________________________
#macro dtm()
For aa As Long=1 To m
For bb As Long=1 To m
If aa=bb Then
det=det *b(aa,bb) 'multiply through the trace of b()
End If
Next bb
Next aa
#endmacro
'The deteminant is found
'___________________________________________________________
#macro invert()
for x2=1 to m
For y2 =1 To m
h(y2)=h2(x2,y2)
Next y2
ivector()
fvector()
For z As Long=1 To m
inverse(z,x2)=solution(z)
Next z
next x2
dtm()
cleanup()
Return sign*det
#endmacro
'**************************** MAIN ***************************
setup()
pivot()
lu()
invert()
'************************** END MAIN *************************
End Function
'Helper subs for examples
Sub matrixmult(m1() As complex,m2() As complex,ans() As complex)
Dim rows As Long=Ubound(m1,1)
Dim columns As Long=Ubound(m2,2)
If Ubound(m1,2)<>UBound(m2,1) Then
Print "Can't do"
Exit Sub
Endif
Redim ans(1 To rows,1 To columns)
Dim rxc As complex
dim null as complex
For r As Long=1 To rows
For c As Long=1 To columns
rxc=null
For k As Long = 1 To Ubound(m1,2)
rxc=rxc+m1(r,k)*m2(k,c)
Next k
ans(r,c)=rxc
Next c
Next r
End Sub
Function rr(byval first As Double, byval last As Double) As Double 'random generator
Function = Rnd * (last - first) + first
End Function
Function decround (byval a As Double,byval b As Integer) As Double 'rounds to b places
Dim y As Double
Dim i As Double
Dim r As Double
y = (Abs(a) - Int(Abs(a))) * (10 ^ b)
i = Int(y)
y = y - i
If y >= .5 Then i = i + 1
i = i / (10 ^ b)
r = Int(Abs(a)) + i
If a < 0 Then r = -r
decround = r
End Function