'Subject: Approximates a real root of integral polynomial f(x)
'         with Lagrange's 1769 continued fractions method.
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DECLARE SUB Evalf (p, q)
'Horner evaluation, let p:= f(q)
DECLARE SUB Coeff ()
'print coefficient vector f[n..0]
CONST a = 0, p0 = 1, p1 = 2, q0 = 3, q1 = 4
CONST e = 100 '                        customize precision e
CONST Fl = 0 '                         verbosity flag
un = LargeInit(q1 + 32, "cnfcroot")
DIM SHARED m(31), n
DIM g AS STRING
CLS : tim# = TIMER: lin = 9

begin:
LOCATE 1, 1
PRINT "   +----------------------------------+"
PRINT "   |  To solve f(x) = 0, input        |"
PRINT "   |  integral coefficients f[n;;0],  |"
PRINT "   |  and the largest integer <= x,   |"
PRINT "   |  f = 0 to exit.                  |"
PRINT "   |                                  |"
PRINT "   +----------------------------------+"
DO: LOCATE 6, 7: INPUT "f"; g
LOOP WHILE INSTR(g, "'")
IF VAL(g) = 0 THEN GOTO eind
LOCATE lin, 1

g = LTRIM$(RTRIM$(g)) + ";"
t = q1: i = 1: k = LEN(g)
DO: t = t + 1
   IF t - q1 > 32 THEN '               maximum rank
      Prints "overflow", 2: GOTO eind
   END IF
   j = INSTR(i, g, ";")
   Readst t, MID$(g, i, j - i) '       read coefficients
   i = j + 1
LOOP UNTIL i > k
n = t - q1 - 1 '                       degree n
FOR i = 0 TO n: m(i) = t - i: NEXT '   initialize f

LOCATE 6, 7: PRINT SPACE$(32); "|"
LOCATE 6, 7: INPUT ; "[x]"; g
Readst a, g '                          read initial a

LOCATE lin, 1: Coeff
IF n = 0 THEN
   Prints "not a polynomial", 2: GOTO begin
END IF

Letf p0, 0: Letf p1, 1 '               initialize recurrences
Letf q0, 1: Letf q1, 0
t = 0
DO: t = t + 1
   Evalf t0, a: i = Isf(t0, 0)
   DO
      Inc a, 1
      Evalf t1, a: j = Isf(t1, 0)
      IF i OR j THEN '                 root found
         IF i THEN Inc a, -1
         IF j THEN Swp t0, t1
         EXIT DO
      END IF
      IF Gets(t0) <> Gets(t1) THEN '   root bracketed
         Inc a, -1: EXIT DO
      END IF
      Swp t0, t1: i = j '              t0:= t1
      IF Getl(a) > 1 THEN '            partial quotients < 32,768
         Prints "search aborted", 2 '  bracketing fails if there
         GOTO begin '                  are 2k roots in [a, a + 1[
      END IF
   LOOP
   Sw = Isf(t0, 0)
   Printn a, "", ",", 0
   '                                   Euclid's algorithm:
   Mult p1, a, t2: Add p0, t2 '        p_t = p_t-2 + a_t * p_t-1
   Mult q1, a, t2: Add q0, t2 '        q_t = q_t-2 + a_t * q_t-1
   Swp p0, p1: Swp q0, q1
   '
   FOR i = 0 TO n - 1
      FOR j = n - 1 TO i STEP -1
         Mult a, m(j + 1), t2 '        Horner substitution y:= x - a
         Add m(j), t2
      NEXT j
   NEXT i
   FOR i = 0 TO (n - 1) \ 2 '          substitute x:= 1 / y
      Swp m(i), m(n - i) '             and multiply by highest power
   NEXT i
   Letf a, 1 '                         reset a
   IF Fl THEN
      Prints "", 1: Coeff '            display all transformed polynomials
   END IF
   i = Getl(p1) + 2 > un '             precision exhausted
   j = Bitl(q1) * .60206 > e
   j = j OR Sw '                       error < 1/ q^2 or root found
LOOP UNTIL i OR j
Prints " [" + LTRIM$(STR$(t)) + "]", 1 ' length of convergent

Printn p1, " p = ", "", 2
Printn q1, " q = ", "", 2
Printr p1, q1, " r = ", e
Prints "", 1
lin = CSRLIN: GOTO begin

eind:
LOCATE lin, 1
g = STR$(CSNG(TIMER - tim#))
Prints "Timer:" + g + " s", 1
CALL Term
SYSTEM

SUB Coeff
   Prints "C: [", 0
   FOR i = n TO 0 STEP -1
      Printn m(i), " ", " ", 0
   NEXT i
   Prints "] d =" + STR$(n), 1
END SUB

SUB Evalf (p, q)
   Dup m(n), p
   FOR i = n - 1 TO 0 STEP -1
      Lmul p, q: Add p, m(i)
   NEXT i
END SUB

