Polynomial Plotter

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Polynomial Plotter

Post by dodicat »

Somebody might find a use for this, but I have my doubts.
I'll maybe try some sort of gui interface, as an exercise, but in the meantime:
edited: gui interface

Code: Select all


'POLY PLOTTER and roots  edit 3rd November
#define WIN_INCLUDEALL
#Include "windows.bi"
#include "fbgfx.bi"
#include "string.bi"
type complex
    as double re,im
end type

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
' USAGE
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"
#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 ********************

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)

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
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
'maxx=maxx+.1*abs(maxx-minx)'extend the x and y ranges a bit for first plot
'minx=minx-.1*abs(maxx-minx)
'maxy=maxy+.1*abs(maxy-miny)
'miny=miny-.1*abs(maxy-miny)
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

Last edited by dodicat on Feb 27, 2012 19:21, edited 1 time in total.
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

I can add complex numbers to Clasp. I have to add vectors first though.
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

rolliebollocks wrote:I can add complex numbers to Clasp. I have to add vectors first though.
Hi again Rollie
I ran your clasp a couple of days ago.
I usually try and keep up with your posts, scattered as they may be among threads.
I suppose Lisp is one of the heavyweights in programming, and is used in many commercial software packages.

I am getting no feedback for my matrix or polynomial stuff, so whether it works on other computers or not, I don't know, but I guess it would.

I thought that for such a large forum there would be somebody could run it and give a Yea or a Nea.

Like your own clasp and library, it takes a bit of messing around to get things up and running, but in my case, it is hardly worth the trouble any more posting maths stuff.
So, in view of that, I'm going to pack in the smoking and draw a breath of fresh air.
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

Part of that file was cut and pasted from a QB document. It's tough to read, and it requires specialist knowledge to appreciate. To a layman it's tough to understand. The user input is really rough. You could use BasicScience's text input, or make a text box yourself.

It would be good to see the coefficients and the polynomial on the screen with the plot and to adjust them dynamically.

You should try the code beutifier or something. Not to sound insulting. But your code is not beautiful.

I need to take a break on programming myself thugh. I'll finish the strings tomorrow.
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

Hi Rollie
I've got this plotter into a gui.
If you follow the instructions, it works ok.
It is not very robust at the moment, I need to fiddle around with the windows gui stuff, it's a bit new to me.
I kinda followed Albert's code to make the gui window, so BLAME ALBERT if all is not well.
edit:
The finished plotter is in the first post.
Last edited by dodicat on Feb 27, 2012 19:25, edited 2 times in total.
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

@Dodicat

Much nicer visually. Still I think the interface is a little awkward. There are too many instructions. People want the input fields to be as obvious as possible. There was a moment when it felt like it was about to crash, but just stalled a while... That should be fixable with a sleep 1 maybe.

What are the 3 different iterations for? You hit space 3 times after hitting plot.

Anyway, definitely headed in the right direction.
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

rolliebollocks wrote:@Dodicat

Much nicer visually. Still I think the interface is a little awkward. There are too many instructions. People want the input fields to be as obvious as possible. There was a moment when it felt like it was about to crash, but just stalled a while... That should be fixable with a sleep 1 maybe.

What are the 3 different iterations for? You hit space 3 times after hitting plot.

Anyway, definitely headed in the right direction.
Hi Rollie
Yes, I've done away with the press key iterations, now do them from scaleup button in the gui.
Also, I've written the roots directly to the gui, it seems a bit better that way.
Sleep doesn't seem to work in a gui loop, a bit like sleep in between screenlock and screenunlock, but it is not needed anyway now.
I've edited the code, just be carefull to enter a number as it should be, no fancy +56.88 or the like.
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Post by dodicat »

Hopefully this is the last version of this.

Code: Select all


'POLY PLOTTER and roots  edit 3rd November
#define WIN_INCLUDEALL
#Include "windows.bi"
#include "fbgfx.bi"
#include "string.bi"
type complex
    as double re,im
end type

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
' USAGE
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"
 #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 ********************

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)

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
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
'maxx=maxx+.1*abs(maxx-minx)'extend the x and y ranges a bit for first plot
'minx=minx-.1*abs(maxx-minx)
'maxy=maxy+.1*abs(maxy-miny)
'miny=miny-.1*abs(maxy-miny)
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
 

 
D.J.Peters
Posts: 8586
Joined: May 28, 2005 3:28
Contact:

Post by D.J.Peters »

good job

Joshy
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york

Post by rolliebollocks »

That turned out really well. GUI looks good, graph looks good. Nice one.
Post Reply