## Polynomial Plotter

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

### Polynomial Plotter

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
else
end if
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=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
if abs(miny)<abs(maxy) then
else
end if
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=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
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
I can add complex numbers to Clasp. I have to add vectors first though.
dodicat
Posts: 7395
Joined: Jan 10, 2006 20:30
Location: Scotland
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
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: 7395
Joined: Jan 10, 2006 20:30
Location: Scotland
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
@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: 7395
Joined: Jan 10, 2006 20:30
Location: Scotland
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: 7395
Joined: Jan 10, 2006 20:30
Location: Scotland
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
else
end if
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=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
if abs(miny)<abs(maxy) then
else
end if
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=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
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: 8443
Joined: May 28, 2005 3:28
Contact:
good job

Joshy
rolliebollocks
Posts: 2655
Joined: Aug 28, 2008 10:54
Location: new york
That turned out really well. GUI looks good, graph looks good. Nice one.