FBComplex: complex number library

Headers, Bindings, Libraries for use with FreeBASIC, Please include example of use to help ensure they are tested and usable.
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

FBComplex: complex number library

Post by jdebord »

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.
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Post by Richard »

Excellent, now when I try to constructively criticise it, all I can come up with is;
Compiling fbcomplex.bas with -w pedantic throws about 6 dozen warnings.
warning 16(0): No explicit BYREF or BYVAL, at parameter .......
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FBComplex: complex number library

Post by dodicat »

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.
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Post by jdebord »

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 :)
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

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!
SARG
Posts: 1756
Joined: May 27, 2005 7:15
Location: FRANCE

Post by SARG »

@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.
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Post by Richard »

Being -w pedantic as ever, Line 73 of testcomp.bas needs a Byval added.
It becomes; SUB TestFunc(Byval Index AS INTEGER)
Well done.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

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.
Sorry SARG, and I was born in 48, probably too young to pick up on it.

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
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Post by jdebord »

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:

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
compiles fine, but when I add either BYVAL or BYREF before Coef() I get an error :

Code: Select all

Illegal specification, at parameter 2 (Coef) of CPoly()
I thought that arrays are passed by reference.
fxm
Moderator
Posts: 12081
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

Post by fxm »

jdebord wrote:I thought that arrays are passed by reference.
Extract of manual about keyword 'Function':
...
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...
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

jdebord 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:

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
compiles fine, but when I add either BYVAL or BYREF before Coef() I get an error :

Code: Select all

Illegal specification, at parameter 2 (Coef) of CPoly()
I thought that arrays are passed by reference.
Hi jdebord
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


 
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

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.

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


jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Post by jdebord »

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:

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
jdebord
Posts: 547
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Post by jdebord »

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 
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

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.


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
 
Post Reply