'Subject: Power basis representation of the fundamental unit
'         (p + qûd) / 2 for the quadratic ring K = Z(ûd).
'         Also given are the continued fraction expansion of the
'         irrational and a decimal approximation of the real square root.
'Author : Sjoerd.J.Schaper
'Date   : 13-03-2005
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST d = 0, p = 1, q = 2, a = 3, b = 4
CONST p0 = 5, p1 = 6, q0 = 7, q1 = 8
t = LargeInit(q1, "fundunit")
DIM g AS STRING, h AS STRING
CLS : tim# = TIMER: lin = 7
CONST dl = 456 '                       decimal length

begin:
LOCATE 1, 1
PRINT "   +---------------------+"
PRINT "   |  input radicand d,  |"
PRINT "   |  0 to exit.         |"
PRINT "   |                     |"
PRINT "   +---------------------+"
LOCATE 4, 7: INPUT ; "d"; di$
IF INSTR(di$, "'") THEN GOTO begin
IF VAL(LEFT$(di$, 5)) < 2 THEN GOTO eind

LOCATE lin, 1
Readst d, di$
Printn d, " d = ", "", 1

Dup d, b: t = 0
Triald p, b, 1024 '                    is D quadratfrei?
WHILE ExQ(b, k, p, t)
   IF k > 1 THEN
      Squ b, p
      Printn p, " square factor ", "", 1
      Prints "", 1: GOTO begin
   END IF
WEND

Isqrt d, b: Dup b, p '                 initial reduced form
iQ = 1: s = Mp2(p, 2) '                store Q_0 and parity(P_0)
IF Mp2(d, 4) = 1 THEN '                D = 1 (mod 4)
   IF s = 0 THEN Inc p, -1
   iQ = 2: s = 1
END IF
Letf q, CLNG(iQ)

Letf p0, 0: Letf p1, 1
Letf q0, 1: Letf q1, 0: t = 0
DO
   Add p, b: t = t + 1
   Divd p, q, a '                      partial quotient a:= (P + sqrt(D)) / Q
   '
   IF t = 1 THEN
      Dup b, a '                       initial a:= sqrt(D), or
      IF Mp2(d, 4) = 1 THEN '           a:= (sqrt(D) - 1) / 2
         Inc a, -1: Rsft a, 1
      END IF
   END IF
   Printn a, "", ",", 0
   '                                   Euclid's algorithm
   Mult p1, a, t0: Add p0, t0 '        p_t = p_t-2 + a_t * p_t-1
   Mult q1, a, t0: Add q0, t0 '        q_t = q_t-2 + a_t * q_t-1
   Swp p0, p1: Swp q0, q1
   '
   Subt p, b: Chs p '                  P':= a * Q - P
   r = Mp2(p, 2) = s '                 is parity(P) = parity(P_0)
   Squ p, t0
   Subt t0, d: Chs t0
   Divd t0, q, t1: Swp q, t1 '         Q':= (D - P'^2) / Q
LOOP UNTIL r AND Isf(q, iQ) '          is Q = Q_0
Prints " [" + LTRIM$(STR$(t)) + "]", 1 ' full period

g = " e = ": r = 0
IF Mp2(d, 4) = 1 THEN  '               power basis representation
   Lsft p1, 1: Add p1, q1 '            x = 2 * p_t + q_t, y = q_t
   r = Mp2(p1, 2) = 1 '                parity(x) = parity(y)
   IF r THEN
      g = g + "("
   ELSE
      Rsft p1, 1: Rsft q1, 1
   END IF
END IF
Printn p1, g, " + ", 0: g = ""
IF NOT Isf(q1, 1) THEN
   Printn q1, "", "", 0: g = " x "
END IF
g = g + "sqrt(": h = ")"
IF r THEN h = h + ") / 2"
Printn d, g, h, 0
Prints "", 1

k = INT(dl * 3.322): t = 0
WHILE Bitl(q1) < k '                  approximate real root(D)
   Squ p1, p: Squ q1, q
   Mult q, d, t0: Add p, t0 '          p_t^2 + q_t^2 * D
   Lmul q1, p1: Lsft q1, 1 '           2p_t * q_t
   Swp p1, p: t = t + 1
WEND
g = " p/q [e^(2^" + LTRIM$(STR$(t)) + ")] = "
Printr p1, q1, g, dl
Prints "", 1
lin = CSRLIN: GOTO begin

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

