'Subject: Lucas pseudoprime test, quick calculation of Fibonacci giants.
'Author : Sjoerd.J.Schaper
'Date   : 05-06-2004
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DECLARE SUB Lucadd (x, y, r)
'Add x and y using small invariant P
DECLARE SUB Lucdbl (x, r)
'Double x using parameter r
DECLARE SUB Lucnum (x, k, N)
'Calculate Lucas number U_k(P, Q) mod N

CONST Q = -1 '                        set parameter Q =+/-1
DIM g AS STRING
DIM SHARED D, y, P, sw
N = 0: D = 1: x = 2: y = 3: k = 4 '   largeint pointers
e = LargeInit(k, "lucasfun")
CLS : tim# = TIMER: lin = 8
IF ABS(Q) <> 1 THEN GOTO eind

begin:
LOCATE 1, 1
PRINT "   +------------------------------------------------------+"
PRINT "   |  Input P, k, (N) to calculate U_k(P,Q) (modulo N)    |"
PRINT "   |  Input k = 'N-e' to perform Lucas pseudoprime-test,  |"
PRINT "   |  P = 0 to exit.                                      |"
PRINT "   |                                                      |"
PRINT "   +------------------------------------------------------+"
DO: LOCATE 5, 7: INPUT ; "P"; g
LOOP WHILE INSTR(g, "'")
IF VAL(g) < 1 OR VAL(g) > MD THEN GOTO eind
P = VAL(g): Readst x, g '             parameter P
INPUT ; " k"; g
fl = INSTR(UCASE$(g), "N") > 0  '     flag to apply primality testing
Readst k, g: Sets k, 1 '              rank abs(k)
INPUT " N"; g
Readst N, g: Sets N, 1 '              modulus abs(N)
g = "": sw = NOT Isf(N, 0) '          switch to apply reduction mod N
IF fl THEN Dup N, k

LOCATE lin, 1: Letf y, 5
Squ x, D: Inc D, -4 * Q '             discriminant D
IF Cmp(D, y) = -1 THEN
   Prints " D < 5", 2: GOTO begin
END IF

IF sw THEN
   e = Kronec(D, N) '                 quadratic residuosity
   IF e THEN
      Euclid D, N, t3 '               D^ -1 (mod N)
      IF fl THEN
	 Inc k, -e
	 IF e = 1 THEN
	    g = "  k = N - 1"
	 ELSE
	    g = "  k = N + 1"
	 END IF
      END IF
   ELSE
      g = " (D/N)= 0"
   END IF
END IF

Printn x, "  P = ", g, 1
IF fl = 0 THEN Printn k, "  k = ", "", 1
IF sw THEN
   Printn N, "  N = ", "", 2
   IF e = 0 THEN Prints "", 1: GOTO begin
END IF

Lucnum x, k, N

IF sw THEN
   Printn x, "U_k ð ", " mod N", 1
   IF fl THEN '                       evaluate primality test
      IF Isf(x, 0) THEN
	 Prints " N is a Lucas pseudoprime", 1
      ELSE
	 Prints " N is composite", 1
      END IF
   END IF
ELSE
   Printn x, "U_k = ", "", 2
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 Lucadd (x, y, r)
   Lmul x, y
   Inc x, -P * r
END SUB

SUB Lucdbl (x, r)
   Squ x, t0: Swp x, t0
   Inc x, -2 * r
END SUB

SUB Lucnum (x, k, N)
   IF Isf(k, 0) THEN
      Letf x, 0: EXIT SUB
   END IF
   q0 = 1: q1 = 1
   Swp y, x: Letf x, 2
   i = Getl(k) - 1
   FOR j = i TO 0 STEP -1 '           L=>R addition chain
      a = Getw(k, j): a0 = MB \ 2
      IF j = i THEN a0 = Hp2(k)
      WHILE a0
	 q0 = q0 * q1
	 IF a AND a0 THEN
	    Lucadd x, y, q0
	    q1 = Q * q0: Lucdbl y, q1
	 ELSE
	    Lucadd y, x, q0
	    Lucdbl x, q0: q1 = q0
	 END IF
	 IF sw THEN
	    Divd x, N, t0
	    Divd y, N, t0
	 END IF
	 a0 = a0 \ 2
      WEND
   NEXT j
   '
   Letf t0, CLNG(P)
   Mult x, t0, t1
   Lsft y, 1: Subt y, t1
   IF sw THEN
      Mult y, D, x
      Modbal x, N
   ELSE '                             U_k = (2V_k+1 - PV_k) / D
      Divd y, D, x
   END IF
END SUB

