10    rem - cubicsolve - a program to find roots to cubic equations
20    rem - with real coeffecients - j.h.stacey - keynes
30    print "entering cubicsolve - 6/11/78 version"
40    print "note this program is not appropriate where roots occur"
50    print "near to the ordinate (x) axis."
60    print
70    print "for your equation ax^3+bx^2+cx+d, input a,b,c,d"
80    input a,b,c,d
90    if a = 0 then 560
100   if d = 0 then 500
110   if a > 0 then 160
120   let a = -a
130   let b = -b
140   let c = -c
150   let d = -d
160   def fna(x) = a*x^3+b*x^2+c*x+d
170   def fnb(x) = 3*a*x^2+2*b*x+c
180   let x1 = 0
190   if b^2 <= 3*a*c then 270
200   let t1 = (-b-sqr(b^2-3*a*c))/(3*a)
210   let t2 = (-b+sqr(b^2-3*a*c))/(3*a)
220   let x1 = (5*t1-3*t2)/2
230   if fna(t2) > 0 then 270
240   let x1 = (5*t2-3*t1)/2
250   if fna(t1) < 0 then 270
260   let x1 = (t1+t2)/2
270   for i = 1 to 50
280    if fnb(x1) = 0 then 430
290    let x = x1-fna(x1)/fnb(x1)
300    if abs(x-x1)/(abs(x)+abs(x1)) < 5e-8 then 430
310    if i = 1 then 370
320    if abs(x-x1) < abs(x1-x2) then 370
330    print "not converging, continue ?";
340    input a$
350    if a$ = "yes" then 370
360    if a$ <> "y" then 730
370    let x2 = x1
380    let x1 = x
390   next i
400   print "convergence criteria not met, ";
410   print "first root set to zero, as implied"
420   let x = 0
430   if x = 0 then 460
440   let z = int(6-log(abs(x))/(log(10)*.999999))
450   let x = sgn(x)*int(abs(x)*10^z+.5)*10^(-z)
460   let e = a
470   let f = x*a+b
480   let g = x*f+c
490   goto 540
500   let x = 0
510   let e = a
520   let f = b
530   let g = c
540   print "solutions are x = ";x;" , ";
550   goto 670
560   if b = 0 then 620
570   print "equation is quadratic, solutions are x = ";
580   let e = b
590   let f = c
600   let g = d
610   goto 670
620   if c = 0 then 650
630   print "equation is linear, soltn is x = ";-d/c
640   goto 730
650   print d;" = 0 is insoluble"
660   goto 730
670   let z = sqr(abs(f^2-4*e*g))
680   if f^2 < 4*e*g then 710
690   print (z-f)/(2*e);" , ";(-f-z)/(2*e)
700   goto 730
710   print -f/(2*e);"+";z/(2*e);"j &"
720   print ,,-f/(2*e);"-";z/(2*e);"j"
730   print "do you have another eqtn. to solve ?";
740   input a$
750   if a$ = "yes" then 80
760   if a$ = "y" then 80
770   print "exiting cubicsolve"
780   end