Back2BASIC

You are here: Home > Issues > Issue #7 > B2B Code Show: Polyplotter

B2B Code Show: Polyplotter

by dodicat

From dodicat:

 

Basically the roots of a real polynomial can be found by iterating,
by QR algorithm, for the eigenvalues of the companion matrix of the
polynomial. The roots are shown with any real roots marked real, of course. The polynomial is sketched, it can be scaled up by a click on the  interface, the roots are drawn onto the sketch.  New coefficients for a different polynomial can be edited at runtime.
'POLY PLOTTER and roots (AUG 2012) 
#define WIN_INCLUDEALL
#Include "windows.bi"
#include "fbgfx.bi"
#include "string.bi"
Type complex
    As Double re,im
End Type
'======================= All procedures =========================
Declare Sub balance(A()As Double,d() As Double)
Declare Sub GET_HESSENBERG(A() As Double)
Declare Sub ITERATE_BY_QR(A() As Double,wr() As Double,wi() As Double)
Declare Sub GET_EIGENVALUES(A() As Double,Eigen() As complex)
Declare Sub makecompanion(polynomium() As Double,mat() As Double)
Declare Sub equalize(m1() As Double,m2() As Double)
Declare Function decround ( a As Double, b As Integer) As Double 
Declare Function _poly(coff() As Double,number As Double)As Double
Declare Function lenint(num As Ulongint) As Integer
Declare Sub getrange(n As Integer)
Declare Sub drawplot(k As Double=1)
Declare Sub plotclick()
Declare Sub scaleupclick
Declare Function seeterms(pol() As Double)As String
'=========  Globals ========================
Dim Shared n As Integer                        'DIMENSION OF POLYNOMIAL
Dim Shared As Double PLOT_GRADE =20000
Dim Shared As Integer xres,yres
Dim Shared As Double minx,maxx,miny,maxy,x
Dim Shared runplot As Integer
Screen 19,32,FB.GFX_ALWAYS_ON_TOP
Screencontrol FB.SET_WINDOW_POS,0,0 
Screeninfo xres,yres
Windowtitle "POLYNOMIAL PLOTTER"
'============== Plotting macros ======================
#macro _window(topleftX,topleftY,bottomrightX,bottomrightY) 
minx=topleftX
maxx=bottomrightX
miny=bottomrightY
maxy=topleftY
#endmacro
#macro _axis(colour)
Line(0,(yres-(miny/(miny-maxy))*yres))-(xres,(yres-(miny/(miny-maxy))*yres)),colour 'x axis
Line(((minx/(minx-maxx))*xres),0)-(((minx/(minx-maxx))*xres),yres),colour 'y axis
Draw String(0,(yres-(miny/(miny-maxy))*yres)),Format(decround(minx,2)),colour
Draw String(xres-8-8*(Len(Format(decround(maxx,2)))),(yres-(miny/(miny-maxy))*yres)),Format(decround(maxx,2)),colour
Draw String (xres/2,0),Format(decround(maxy,2)),colour
Draw String (xres/2,yres-16),Format(decround(miny,2)),colour
#endmacro
#macro grid(colour)
Scope
    Dim As Integer l
    Dim As Double grader
    If Abs(minx)<Abs(maxx) Then
        grader=minx
    Else
        grader=maxx
    End If
    If grader<=18446744073709551615 Then
        l=Int(Abs(grader))'min
        l=lenint(l)
        For z As Double=0 To minx Step Sgn(minx)*(10^(l-1))
            Line((((minx-z)/(minx-maxx))*xres),0)-((((minx-z)/(minx-maxx))*xres),yres),colour 'y axis 
        Next z
        l=Int(Abs(grader))'max
        l=lenint(l)
        For z As Double=0 To maxx Step Sgn(maxx)*(10^(l-1))
            Line((((minx-z)/(minx-maxx))*xres),0)-((((minx-z)/(minx-maxx))*xres),yres),colour 'y axis 
        Next z
    End If'grader
    If Abs(miny)<Abs(maxy) Then
        grader=miny
    Else
        grader=maxy
    End If
    If grader<=18446744073709551615 Then
        l=Int(Abs(grader))'miny
        l=lenint(l)
        For z As Double=0 To miny Step Sgn(miny)*(10^(l-1))
            Line(0,(yres-((miny-z)/(miny-maxy))*yres))-(xres,(yres-((miny-z)/(miny-maxy))*yres)),colour '    
        Next z
        l=Int(Abs(grader))'max
        l=lenint(l)
        For z As Double=0 To maxy Step Sgn(maxy)*(10^(l-1))
            Line(0,(yres-((miny-z)/(miny-maxy))*yres))-(xres,(yres-((miny-z)/(miny-maxy))*yres)),colour '    
        Next z
    End If
End Scope
#endmacro
#macro sketch(_function,colour)
For x =minx To maxx Step (maxx-minx)/PLOT_GRADE
    Dim As Double x1=Cdbl(xres)*(x-minx)/(maxx-minx)
    Dim As Double y1=Cdbl(yres)*(_function-maxy)/(miny-maxy)
    Pset(x1,y1),colour
Next x
#endmacro
#macro mark(x1,y1,colour)
Dim As Double xx1= Cdbl(xres)*(x1-minx)/(maxx-minx)
Dim As Double yy1=Cdbl(yres)*(y1-maxy)/(miny-maxy)
Circle (xx1,yy1),5,colour,,,,f
#endmacro
' ***************** END PLOTTING MACROS ********************
'=============== Set up Win Gui =============================
Paint (0,0),Rgb(200,200,200)
Dim As MSG msg
Dim Shared As HWND hWnd,hwnd2,Poly_input1,Poly_input2,plot,roots_output,note,scaleup,scaledown
Dim Shared As String outtext,errorflag,rootstring
errorflag=""
Redim Shared As Double poly(0),copypoly(0)
Dim Shared As String s3
Dim dx As Integer=250
'*******************************************************
hWnd=CreateWindowEx( 0, "#32770", "Coefficients", WS_OVERLAPPEDWINDOW Or WS_VISIBLE, 200+dx, 200, 500+dx, 600, 0, 0, 0, 0 )
Dim s As String= "Enter Coefficients, highest power first."+Chr(13)+Chr(10)
Dim s2 As String="When all entered, click PLOT"
s=s+s2
roots_output = CreateWindowEx( 0, "EDIT", "", ws_border Or WS_VISIBLE Or WS_CHILD Or WS_HSCROLL Or WS_VSCROLL Or ES_AUTOHSCROLL Or ES_AUTOVSCROLL Or ES_MULTILINE Or ES_READONLY, 220, 80, 450, 400, hWnd, 0, 0, 0 )
note = CreateWindowEx( 0, "EDIT",s, ws_border Or WS_VISIBLE Or WS_CHILD Or WS_HSCROLL Or WS_VSCROLL Or ES_AUTOHSCROLL Or ES_AUTOVSCROLL Or ES_MULTILINE Or ES_READONLY, 110, 20, 300, 60, hWnd, 0, 0, 0 )
Poly_input2 = CreateWindowEx( 0, "EDIT","", ws_border Or WS_VISIBLE Or WS_CHILD Or WS_HSCROLL Or WS_VSCROLL Or ES_AUTOHSCROLL Or ES_AUTOVSCROLL Or ES_MULTILINE, 10, 80, 200, 400, hWnd, 0, 0, 0 )
plot= CreateWindowEx( 0, "BUTTON", "Plot", WS_VISIBLE Or WS_CHILD, 20, 20 , 70, 30, hWnd, 0, 0, 0 )
scaleup=CreateWindowEx( 0, "BUTTON", "Scaleup", WS_VISIBLE Or WS_CHILD, 20, 20+500 , 70, 30, hWnd, 0, 0, 0 )
s3="INSTRUCTIONS for use:"+Chr(13)+Chr(10)
s3=s3+"Print coefficients carefully and realistically"+Chr(13)+Chr(10)
s3=s3+"Example -3, not -03 or not -0 for 0 or not 0.8 for .8"+Chr(13)+Chr(10)
s3=s3+"No un-needed + signs"+Chr(13)+Chr(10)
s3=s3+"E.g. +.78"+Chr(13)+Chr(10)
s3=s3+"No numbers in form 3.5e-5, you must write .000035 in this case"+Chr(13)+Chr(10)
s3=s3+"Press enter at each new entry, "
s3=s3+"including the last entry"+Chr(13)+Chr(10)
s3=s3+"Edit previous entries if required"+Chr(13)+Chr(10)
s3=s3+"No need to press enter in this case"+Chr(13)+Chr(10)
s3=s3+"First coefficient (highest power) must not be 0"
setWindowText(roots_output,s3)
Dim As String starttext
starttext="16"+Chr(13)+Chr(10)+"0"+Chr(13)+Chr(10)+"-20"+Chr(13)+Chr(10)
starttext=starttext+"0"+Chr(13)+Chr(10)+"5"+Chr(13)+Chr(10)+"0"+Chr(13)+Chr(10)
setwindowtext(poly_input2,starttext)
Draw String (50,50),"Start example is Chebyshev polynomial degree 5 ",Rgb(0,0,200)
Draw String (50,70),"Please read the instructions, then click plot to start",Rgb(0,0,200)
'=======================  Working Loop ==========================
start:
While GetMessage( @msg, 0, 0, 0 )
    TranslateMessage( @msg )
    DispatchMessage( @msg )
    Select Case msg.hwnd
    Case hWnd
        Select Case msg.message
        Case 273
            End
        End Select
        '__________________  
    Case plot
        Select Case msg.message  
        Case WM_LBUTTONDOWN
            plotclick()
            Goto go
        End Select
        '______________________   
    Case scaleup
        Select Case msg.message  
        Case WM_LBUTTONDOWN
            If runplot>=1 Then
                scaleupclick()
            End If
        End Select 
        '__________________________
        
        
    End Select
Wend
'==================== PROCEDURES ============================
Sub scaleupclick
    Static k As Double
    k=k+.001
    drawplot(1+k)
End Sub
Sub plotclick()
    #macro remove(s,char)
    Scope
        Dim temp As String
        Dim As Integer asci=Asc(char)
        For i As Long =0 To Len(s)-1 
            If s[i]<>asci Then temp=temp+Chr$(s[i])
        Next i
        s= temp
    End Scope
    #endmacro
    #macro Replace(s,char,newchar)
    Scope
        Dim temp As String=String(Len(s),newchar)
        Dim As Integer asci=Asc(char)
        For i As Long =0 To Len(s)-1 
            If s[i]<>asci Then temp[i]=s[i]
        Next i
        s= temp
    End Scope
    #endmacro
    
    Dim charcount As Integer
    charcount = GetWindowTextLength(Poly_input2)
    outtext = String(charcount," ")
    GetWindowText(Poly_input2,outtext,charcount)
    remove(outtext,Chr(0))
    outtext=Rtrim(outtext,Chr(10))
    outtext=Rtrim(outtext,Chr(13))
    outtext=Rtrim(outtext,Chr(10))
    replace(outtext,Chr(10),"*")
    replace(outtext,Chr(13),"*")
    outtext=outtext+Chr(0)
    Redim poly(charcount)
    Redim copypoly(charcount)
    Dim count As Integer
    For z As Integer=0 To Ubound(poly)
        count=count+1
        outtext=Ltrim(outtext,"+") 'just incase
        poly(z)=Val(outtext)
        copypoly(z)=poly(z)
        outtext=Ltrim(outtext,Format(poly(z)))
        outtext=Ltrim(outtext,"*")
        If outtext=Chr(0) Then
            outtext=Rtrim(outtext,Chr(0))
            Exit For
        End If
    Next z
    Redim Preserve poly(count-1)
    Redim Preserve copypoly(count-1)
    If poly(0)=0 Then
        messagebox(NULL,"The first coefficient should not be 0","ERROR",MB_OK)
        errorflag="ERROR"+Chr(13)+Chr(10)
        errorflag=errorflag+"First coefficient"+Chr(13)+Chr(10)
        errorflag=errorflag+"can't be 0"
    Else
        errorflag=""
    End If
    
    'REVERSE POLY
    Dim As Long lb,ub:lb=Lbound(poly):ub=Ubound(poly)
    For z As Integer=Lb To Int((lb+Ub)/2)
        Swap poly(z),poly(ub+lb-z)
        Swap copypoly(z),copypoly(ub+lb-z)
    Next z
    
    outtext=s3+Chr(13)+Chr(10)
    outtext=outtext+errorflag
    
    setWindowText(roots_output,outtext)
End Sub
go:
n=Ubound(poly)
If n<=0 Then errorflag="no"
Redim companion(0,0) As Double
Redim Shared As complex roots(0)
If errorflag="" Then 
    Dim As Double val1=poly(n)
    For z As Integer=n To 0 Step -1 'prepare to form companion matrix
        poly(z)=poly(z)/val1          'leading coefficient=1
    Next z
    Print
    
    makecompanion(poly(),companion())
    
    GET_EIGENVALUES(companion(),roots())
    getrange(n)
    
    drawplot(.1)
    
End If 'errorflag=""
Goto start
Sub getrange(n As Integer)
    ' get plotting parameters
    maxx=-1e50
    minx=1e50
    maxy=-1e50
    miny=1e50
    For z As Integer=1 To n
        If maxx<roots(z).re Then maxx=roots(z).re 'get x range
        If minx>roots(z).re Then minx=roots(z).re
    Next z
    If maxx=0 And minx=0 Then
        For z As Integer=1 To n
            If maxx<roots(z).im Then maxx=roots(z).im 'if no real roots, get x range on imaginary roots
            If minx>roots(z).im Then minx=roots(z).im
        Next z
    End If
    If maxx=0 And minx=0 Then 'if roots are all totally zero
        Print "DEFAULT X LIMITS"
        maxx=1
        minx=-1
    End If
    
    If minx=maxx Then
        minx=minx-.1*Abs(minx)
        maxx=maxx+.1*Abs(maxx)
    End If
    Dim As Double polyval
    For z As Double=minx To maxx Step (maxx-minx)/1000 'get y range
        polyval=_poly(copypoly(),z)
        If maxy< polyval Then maxy=polyval
        If miny>polyval Then miny=polyval
    Next z
    
    For z As Double=1 To n
        If maxy< roots(z).im Then maxy=roots(z).im 'extend to cover imaginary roots
        If miny>roots(z).im Then miny=roots(z).im
    Next z
    
End Sub
Sub drawplot(k As Double=1)
    
    Dim As Double xhalf,yhalf,range
    xhalf=Abs(maxx-minx)/2
    yhalf=Abs(maxy-miny)/2
    range=Abs(maxx-minx)
    If range<1 Then k=10*k
    maxx=(maxx+xhalf*k)
    minx=(minx-xhalf*k)
    maxy=(maxy+yhalf*k)
    miny=(miny-yhalf*k)
    Cls
    Paint (0,0),Rgb(200,200,200)
    
    runplot=runplot+1
    _window(minx,maxy,maxx,miny)
    
    grid(Rgb(180,180,180))'draw grid
    
    _axis(Rgb(200,0,0))   'draw axis
    Dim As String cplex,polystring
    For z As Integer=n To 0 Step -1  'write coefficients
        polystring=polystring+Format(copypoly(z))+","
    Next z
    Dim As Uinteger col  'write roots
    
    
    Draw String(.01*xres,16),"COEFFICIENTS "+"("+polystring+")"+" (degree "+Str(n)+")"
    
    Draw String(.3*xres,32), "ROOTS marked (blue = real) (white = imaginary)"
    rootstring="ROOTS:"+Chr(13)+Chr(10)
    For x As Integer=1 To n         'seperate real from imaginary
        cplex=" , "+Str(roots(x).im)+" i"
        col=Rgb(255,255,255)
        If  Val(Str(roots(x).im))=0  Then 
            cplex=" (real root)"
            col=Rgb(0,0,200)
        End If
        rootstring=rootstring+Str(roots(x).re)+cplex+Chr(13)+Chr(10)
        setWindowText(roots_output,rootstring) 
    Next x
    setWindowText(note,seeterms(copypoly()))
    Sketch(_poly(copypoly(),x),Rgb(0,200,0))   'draw the curve
    
    For z As Double=1 To n                     'mark the roots
        mark(roots(z).re,roots(z).im,Rgb(255*Abs(Sgn(roots(z).im)),255*Abs(Sgn(roots(z).im)),255))
    Next z
    
End Sub
'START PROCEDURES FOR EIGENVALUES
#macro EXC(m)
d(m) = 1# * j
If j <> m Then
    For ii As Integer = 1 To k
        f = A(ii, j)
        A(ii, j) = A(ii, m)
        A(ii, m) = f
    Next ii
    For ii As Integer = l To n
        f = A(j, ii)
        A(j, ii) = A(m, ii)
        A(m, ii) = f
    Next ii
End If
#endmacro
Sub balance(A() As Double,d() As Double)
    Dim n As Integer = Ubound(A) 
    Dim As Double b = 2    
    Dim As Integer low,ihi,noconv
    Dim As Integer i,j,k,l,m        
    Dim As Double b2,c,f,g,r,s  
    
    
    Dim As Double zero=0
    b2 = b * b
    l = 1
    k = n
    
    lab1: For j = k To 1 Step -1
    r = ZERO
    For i = 1 To j - 1
        r = r + Abs(A(j, i))
    Next i
    For i = j + 1 To k
        r = r + Abs(A(j, i))
    Next i
    If (r = 0) Then
        m = k
        Exc(k)
        k = k - 1
        Goto lab1
    End If
Next j
lab2: For j = l To k
c = ZERO
For i = l To j - 1
    c = c + Abs(A(i, j))
Next i
For i = j + 1 To k
    c = c + Abs(A(i, j))
Next i
If (c = 0) Then
    m = l: Exc(l)
    l = l + 1
    Goto lab2
End If
Next j
low = l
ihi = k
For i = 1 To k
    d(i) = 1#
Next i
lab3:  noconv = 0
For i = l To k
    c = ZERO
    r = c
    For j = l To i - 1
        c = c + Abs(A(j, i))
        r = r + Abs(A(i, j))
    Next j
    For j = i + 1 To k
        c = c + Abs(A(j, i))
        r = r + Abs(A(i, j))
    Next j
    g = r / b
    f = 1#
    s = c + r
    lab4: If c < g Then
    f = f * b
    c = c * b2
    Goto lab4
End If
g = r * b
lab5: If c >= g Then
f = f / b
c = c / b2
Goto lab5
End If
If (c + r) / f < .95 * s Then
    g = 1# / f
    d(i) = d (i) * f
    noconv = 1
    For j = l To n
        A(i, j) = A(i, j) * g
    Next j
    For j = 1 To k
        A(j, i) = A(j, i) * f
    Next j
End If
Next i
If noconv = 1 Then Goto lab3
End Sub
Sub GET_HESSENBERG(A() As Double)
    
    Dim n As Integer=Ubound(a)
    Dim As Integer m,j,i
    Dim As Double y,x,ZERO=0
    If n > 2 Then
        For m = 2 To n - 1
            x = ZERO
            i = m
            For j = m To n
                If Abs(A(j, m - 1)) > Abs(x) Then
                    x = A(j, m - 1)
                    i = j
                End If 'IF Abs
            Next j 'FOR j= m TO n
            If i <> m Then
                For j = m - 1 To n
                    y = A(i, j)
                    A(i, j) = A(m, j)
                    A(m, j) = y
                Next j
                For j = 1 To n
                    y = A(j, i)
                    A(j, i) = A(j, m)
                    A(j, m) = y
                Next j
            End If 'IF i <> m
            If x <> ZERO Then
                For i = m + 1 To n
                    y = A(i, m - 1)
                    If y <> ZERO Then
                        y = y / x
                        A(i, m - 1) = y
                        For j = m To n
                            A(i, j) = A(i, j) - y * A(m, j)
                        Next j
                        For j = 1 To n
                            A(j, m) = A(j, m) + y * A(j, i)
                        Next j
                    End If  'IF y
                Next i  'FOR i
            End If  'IF x
        Next m  'FOR m
    End If 'if n>2
End Sub
Sub ITERATE_BY_QR(A() As Double,wr() As Double,wi() As Double)
    Dim As Double y,r,u,v,mmin
    Dim As String res
    Dim sign As Double
    Dim As Double ZERO=0
    Dim n As Integer=Ubound(A)
    Dim As Double anorm,itsmx,t,its,s,x,p,q,nn,z,w
    Dim As Integer i,j,l,m,k
    Dim As Double aa,bb
    #macro _SIGN(aa,bb)
    If bb < 0 Then
        Sign = -Abs(aa)
    Else
        Sign = Abs(aa)
    End If
    #endmacro
    itsmx = 50
    
    anorm = Abs(A(1, 1))
    For i = 2 To n
        For j = i - 1 To n
            anorm = anorm + Abs(A(i, j))
        Next j
    Next i
    nn = n
    t = ZERO
    
    label4: its = 0
    label2: For l = nn To 2 Step -1
    s = Abs(A(l - 1, l - 1)) + Abs(A(l, l))
    If (s = 0!) Then s = anorm
    If ((Abs(A(l, l - 1)) + s) = s) Then Goto label3
Next l
l = 1
label3: x = A(nn, nn)
If (l = nn) Then
    wr(nn) = x + t
    wi(nn) = ZERO
    nn = nn - 1
Else
    y = A(nn - 1, nn - 1)
    w = A(nn, nn - 1) * A(nn - 1, nn)
    If (l = nn - 1) Then
        p = .5# * (y - x)
        q = p * p + w
        z = Sqr(Abs(q))
        x = x + t
        If q >= ZERO Then
            aa = z: bb = p:  _SIGN(aa,bb)
            z = p + Sign
            
            wr(nn) = x + z
            wr(nn - 1) = wr(nn)
            If z <> ZERO Then wr(nn) = x - w / z
            
            
            wi(nn) = ZERO
            wi(nn - 1) = ZERO
        Else
            wr(nn) = x + p
            wr(nn - 1) = wr(nn)
            wi(nn) = z
            wi(nn - 1) = -z
        End If
        nn = nn - 2
    Else
        If its = itsmx Then
            Print " Pause"
            Print " Too many iterations!"
            Input res
        End If
        If (its = 10) Or (its = 20) Then
            t = t + x
            For i = 1 To nn
                A(i, i) = A(i, i) - x
            Next i
            s = Abs(A(nn, nn - 1)) + Abs(A(nn - 1, nn - 2))
            x = .75# * s
            y = x
            w = -.4375# * s * s
        End If
        its = its + 1
        For m = nn - 2 To 1 Step -1
            z = A(m, m)
            r = x - z
            s = y - z
            p = (r * s - w) / A(m + 1, m) + A(m, m + 1)
            q = A(m + 1, m + 1) - z - r - s
            r = A(m + 2, m + 1)
            s = Abs(p) + Abs(q) + Abs(r)
            p = p / s
            q = q / s
            r = r / s
            If m = 1 Then Goto label1
            u = Abs(A(m, m - 1)) * (Abs(q) + Abs(r))
            v = Abs(p) * (Abs(A(m - 1, m - 1)) + Abs(z) + Abs(A(m + 1, m + 1)))
            If u + v = v Then Goto label1
        Next m
        label1:  For i = m + 2 To nn
        A(i, i - 2) = ZERO
        If i <> (m + 2) Then A(i, i - 3) = ZERO
    Next i
    For k = m To nn - 1
        If k <> m Then
            p = A(k, k - 1)
            q = A(k + 1, k - 1)
            r = ZERO
            If k <> (nn - 1) Then r = A(k + 2, k - 1)
            x = Abs(p) + Abs(q) + Abs(r)
            If x <> ZERO Then
                p = p / x
                q = q / x
                r = r / x
            End If
        End If
        aa = Sqr(p * p + q * q + r * r): bb = p:  _SIGN(aa,bb)'GOSUB 2900
        s = Sign
        If s <> ZERO Then
            If k = m Then
                If l <> m Then A(k, k - 1) = -A(k, k - 1)
            Else
                A(k, k - 1) = -s * x
            End If
            p = p + s
            x = p / s
            y = q / s
            z = r / s
            q = q / p
            r = r / p
            For j = k To nn
                p = A(k, j) + q * A(k + 1, j)
                If k <> (nn - 1) Then
                    p = p + r * A(k + 2, j)
                    A(k + 2, j) = A(k + 2, j) - p * z
                End If
                A(k + 1, j) = A(k + 1, j) - p * y
                A(k, j) = A(k, j) - p * x
            Next j
            If nn < k + 3 Then
                mmin = nn
            Else
                mmin = k + 3
            End If
            For i = 1 To mmin
                p = x * A(i, k) + y * A(i, k + 1)
                If k <> (nn - 1) Then
                    p = p + z * A(i, k + 2)
                    A(i, k + 2) = A(i, k + 2) - p * r
                End If
                A(i, k + 1) = A(i, k + 1) - p * q
                A(i, k) = A(i, k) - p
            Next i
        End If
    Next k
    Goto label2:
End If
End If
If nn >= 1 Then Goto label4:
End Sub
Sub GET_EIGENVALUES(MAT() As Double,Eig() As complex)
    
    Dim n As Integer=Ubound(MAT)
    Redim Eig(n)
    Dim As Double A(n,n)
    equalize(A(),MAT())
    Dim As Double wr(n),wi(n),d(n)
    balance(a(),d())
    GET_HESSENBERG(A())
    ITERATE_BY_QR(A(),wr(),wi())
    For i As Integer = 1 To n
        Eig(i).re= wr(i):Eig(i).im= wi(i)
    Next i
End Sub
Sub makecompanion(polynomium() As Double,mat() As Double) 
    Dim As Integer di=Ubound(polynomium)
    Dim As Double one,_null=0
    one=1:
    Redim mat (1 To di,1 To di) As Double
    
    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
Sub equalize(m1() As Double,m2() As Double)
    For x As Integer=1 To Ubound(m1)
        For y As Integer=1 To Ubound(m1) 
            m1(x,y)=m2(x,y) 
        Next y
    Next x
End Sub
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
Function lenint(num As Ulongint) As Integer
    Dim As Double x=1
    Dim As Integer l=1
    Do 
        x=x*10   
        Select Case num>=x
        Case -1
            l=l+1
        Case 0
            Exit Select
        End Select
    Loop Until x>=18446744073709551615
    Return l
End Function
Function _poly(coff() As Double,number As Double)As Double
    Dim count As Integer                'evaluates the polynomial
    Dim pol As Double
    Dim deg As Integer=Ubound(coff)
    pol = 0
    For count = 1 To DEG + 1
        pol = pol + coff(count-1) * ((number) ^ (count - 1))
    Next count
    _poly = pol
End Function
Function seeterms(pol() As Double)As String
    Redim As Double coff(1 To Ubound(pol)+1)
    For z As Integer=0 To Ubound(pol)'+1
        coff(z+1)=pol(z)
    Next z
    Dim result As String
    result=result+"P(x) = "
    Dim DEG As Integer=Ubound(coff)-1
    Dim sgnstr As String
    For j As Integer=DEG+1 To 1 Step -1
        
        If coff(j) <> 0 Then
            If Sgn(coff(j)) = 1 Then 
                sgnstr = "+"
                result=result+"+"
            End If
            
            If Sgn(coff(j)) = -1 Then 
                sgnstr = ""
                result=result+""
            End If
            If j - 1 = 0 Then 
                result=result+Format(coff(j))
            End If
            
            If Abs(coff(j)) <> 1 Then
                If j - 1 = 1 Then 
                    result=result+Format(coff(j))+"x "
                End If
                
                
            End If
            If Abs(coff(j)) <> 1 Then
                If j - 1 > 1 Then 
                    result=result+Format(coff(j))+"x^"+(Format(j - 1))+" "
                End If
                
            End If
            If coff(j) = 1 Then
                If j - 1 > 1 Then 
                    result=result+"x^"+(Format(j - 1))+" "
                End If
                
                
                If j - 1 = 1 Then 
                    result=result+"x "
                End If
            End If
            
            
            If coff(j) = -1 Then
                If j - 1 = 1 Then 
                    result=result+"-x "
                End If
                
                If j - 1 > 1 Then 
                    result=result+"-x^"+(Format(j - 1))
                End If
            End If
        End If
    Next j
    Return result
End Function
 

 

 

| top |