'Subject: Pseudoprime test using a k-th order linear recurrence
'         having monic f(x) as characteristic polynomial.
'         The sequence is periodic modulo a prime, and is computed
'         by powering modulo N the companion matrix for f(x).
'         Then if N is prime Tr(A^N) = Tr(A) and Det(A^N) = Det(A).
'         "s" Attached to N will print a k-segment of the recursive
'         sequence starting with the n-th term, and ratio Un+1 / Un.
'
'Author : Sjoerd.J.Schaper
'Date   : 21-07-2005
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DEFINT C
DECLARE SUB Pwrmtx (A, p)
'A:= A ^ p
DECLARE SUB Dupmtx (B, A)
'B:= A
DECLARE SUB Mulmtx (A, B)
'A:= A * B
DECLARE SUB Trace (q, A)
'q:= sum of the roots (mod N)
DECLARE SUB Det (q, A)
'q:= product of roots (mod N)
DECLARE SUB Prntmtx (A)
'print matrix A

CONST Xk = 9 '                         maximal degree - 1
CONST Le = 28 '                        decimal fraction length
DIM SHARED M(2, Xk, Xk), k
DIM SHARED N, B, C, red
DIM g AS STRING, c0 AS STRING, ck AS STRING
N = 0: q = 1: i = q '                  largeint pointers
A = 0: B = 1: C = 2 '                  matrix pointers
FOR t = A TO C
   FOR r = 0 TO Xk
      FOR s = 0 TO Xk
	 i = i + 1: M(t, r, s) = i '   initialize indices
      NEXT s
   NEXT r
NEXT t
r = LargeInit(i, "powrmtrx")
CLS : tim# = TIMER: lin = 9

begin:
LOCATE 1, 1
PRINT "   +-------------------------------+"
PRINT "   |  Input integral coefficients  |"
PRINT "   |  f[n-1;;0] of a monic poly,   |"
PRINT "   |  and number N to test,        |"
PRINT "   |  f = 0 to exit.               |"
PRINT "   |                               |"
PRINT "   +-------------------------------+"
DO: LOCATE 6, 7: INPUT "f"; g
LOOP WHILE INSTR(g, "'")
IF g = "" OR g = "0" THEN GOTO eind
LOCATE lin, 1

g = LTRIM$(RTRIM$(g)) + ";"
k = -1: i = 1: r = LEN(g)
DO: k = k + 1 '                        read coefficients
   IF k > Xk THEN
      Prints "overflow", 2
      lin = CSRLIN: GOTO eind
   END IF
   j = INSTR(i, g, ";")
   c0 = MID$(g, i, j - i)
   Readst M(A, 0, k), c0
   IF k = 0 THEN ck = c0
   i = j + 1
LOOP UNTIL i > r

LOCATE 6, 7: PRINT SPACE$(29); "|"
LOCATE 6, 7: INPUT ; "N"; g '          modulus or exponent
IF VAL(g) < 2 THEN GOTO begin
LOCATE lin, 1: red = -1
r = INSTR(LCASE$(g), "s")
IF r THEN
   MID$(g, r, 1) = " ": red = 0 '      disable reduction
   IF VAL(g) > MD THEN
      Prints " overflow", 2: GOTO begin
   END IF
END IF
Readst N, g

Prints "C: [ 1 ", 0
FOR s = 0 TO k
   Printn M(A, 0, s), " ", " ", 0
NEXT s
Prints "] d =" + STR$(k + 1), 1
Printn N, " N = ", "", 1 - red

FOR s = 0 TO k
   Chs M(A, 0, s)
   SWAP M(A, 0, s), M(A, k, k - s) '   build companion matrix
NEXT s
FOR r = 0 TO k - 1
   FOR s = 0 TO k
      IF s - r = 1 THEN
	 Letf M(A, r, s), 1
      ELSE
	 Letf M(A, r, s), 0
      END IF
   NEXT s
NEXT r

Pwrmtx A, N
'Prntmtx A

IF red THEN '                          evaluate primality test
   Trace q, A
   Printn q, " Tr = ", " mod N", 1
   Readst t0, ck: Chs t0
   Modbal t0, N
   sw = Cmp(q, t0) = 0
   '
   Det q, A
   Printn q, " Det = ", " mod N", 1
   Readst t0, c0
   IF (k + 1) AND 1 THEN Chs t0
   Modbal t0, N
   IF sw AND (Cmp(q, t0) = 0) THEN
      Prints "N is pseudoprime relative to f(x)", 1
   ELSE
      Prints "N is composite", 1
   END IF
ELSE
   '
   i = M(A, 0, k)
   Printn i, " S: ", " ", 0
   FOR r = 1 TO k
      Printn M(A, r, k), " ", " ", 0 ' sequence
   NEXT r: Prints "", 1
   IF k > 0 AND NOT Isf(i, 0) THEN '   approximate root
      Printr M(A, 1, k), i, " Un+1/Un = ", Le
   END IF
END IF

Prints "", 1
lin = CSRLIN
GOTO begin

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

SUB Det (q, A)
   Letf q, 1: e = 1
   FOR s = 0 TO k - 1
      FOR r = s TO k
	 Gcd M(A, r, s), N, t3 '       find invertible pivot
	 sw = Isf(t3, 1)
	 IF sw THEN EXIT FOR
      NEXT r
      IF sw = 0 THEN Letf q, 0: EXIT SUB
      '
      IF r > s THEN
	 FOR t = s TO k '              exchange rows,
	    SWAP M(A, s, t), M(A, r, t)
	 NEXT t: e = -e '              change permutation sign
      END IF
      Modmlt q, M(A, s, s), N
      '
      Euclid M(A, s, s), N, t3 '       invert pivot
      FOR r = s + 1 TO k '             Gaussian elimination
	 Mult M(A, s, s), M(A, r, s), t3
	 Divd t3, N, t0
	 FOR t = s + 1 TO k
	    Mult t3, M(A, s, t), t0
	    Divd t0, N, t1
	    Subt M(A, r, t), t0
	 NEXT t
      NEXT r
   NEXT s
   Lmul q, M(A, k, k) '                Det(A)
   IF e = -1 THEN Chs q
   Modbal q, N
END SUB

SUB Dupmtx (B, A)
   FOR r = 0 TO k
      FOR s = 0 TO k
	 Dup M(A, r, s), M(B, r, s)
      NEXT s
   NEXT r
END SUB

SUB Mulmtx (A, B)
   FOR t = 0 TO k
      FOR r = 0 TO k
	 Letf M(C, r, t), 0
	 FOR s = 0 TO k
	    Mult M(A, r, s), M(B, s, t), t0
	    Add M(C, r, t), t0
	 NEXT s
	 IF red THEN Divd M(C, r, t), N, t0
      NEXT r
   NEXT t
   SWAP A, C
END SUB

SUB Prntmtx (A)
   FOR r = 0 TO k
      FOR s = 0 TO k
	 Printn M(A, r, s), " ", " ", 0
      NEXT s: Prints "", 1
   NEXT r: Prints "", 1
END SUB

SUB Pwrmtx (A, p)
   i = Getl(p) - 1
   FOR j = i TO 0 STEP -1 '            L=>R binary exponentiation
      e = Getw(p, j): e0 = MB \ 2
      IF j = i THEN
	 Dupmtx B, A: e0 = Hp2(p) \ 2
      END IF
      WHILE e0
	 Mulmtx B, B
	 IF e AND e0 THEN Mulmtx B, A
	 e0 = e0 \ 2
      WEND
   NEXT j
   SWAP A, B
END SUB

SUB Trace (q, A)
   Letf q, 0
   FOR r = 0 TO k
      Add q, M(A, r, r) '              trace(A)
   NEXT r
   Modbal q, N
END SUB

