rem roots of polynomials - half interval search
cls
print
print "roots of polynomials - half interval search"
print "please enter the equation in 'function f(x)'"
print "to end search enter 'q' for either interval"

dim d(3)
lim1=.000005
lim2=.00001

[interval]
print
input "enter lower interval (q to quit): ";a$
if a$="q" then end
input "enter upper interval (q to quit): ";b$
if b$="q" then end
a=val(a$)
b=val(b$)
print
if a=b or a>b then
 gosub interval
else
 gosub [calc]
end if
goto [interval]

rem ***********************************************************************
rem *****  enter function below, example f=2*x^3-1 *************************
rem ***********************************************************************
function f(x)
 f=4*x^4-2.5*x^2-x+.5
rem f=x^2
end function
rem ***********************************************************************

[calc]
a1=sgn(f(a))
b1=sgn(f(b))
rem is either a root?
if a1=0 then
 y=a
 gosub [print]
end if
if b1=0 then
 y=b
 gosub [print]
end if
if a1*b1<0 then gosub [found]

rem loop to search 3000 times for sign change or zero on interval
for i=1 to 3000
 x=a+rnd(2)*(b-a)
 x2=f(x)
 x1=sgn(f(x))
 if x1=0 then
  y=x
  gosub [print]
 end if
 if ((x2<lim2) and (x2>-1*lim2)) then
   print
   print "possible root at ";x
   goto [interval]
 end if
 if a1*x1<0 then
  b=x
  gosub [found]
 end if
next i
print "no roots found in interval..."
return

[found]
rem store positive point in d3
rem store negative point in d1
d(2+a1)=a
d(2-a1)=b
[loop]
y=(d(1)+d(3))/2
y1=sgn(f(y))
if y1=0 then gosub [print]
d(2+y1)=y
if abs(d(1)-d(3))/abs(d(1)+abs(d(3)))<lim1 then gosub [print]
goto [loop]
return

[print]
print
print "root= ";y
print
goto [interval]
return

function sgn(num)
 if num>0 then sgn=1
 if num<0 then sgn=-1
 if num=0 then sgn=0
end function