'Subject: Least primitive root of N, trial division up to 2^ 20.
'         Integer N has a primitive root if N = 2, 4, p^ k or 2p^ k.
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST a = 0, b = 1, N = 2, p = 3, q = 4, R = 5
s = LargeInit(R, "primroot")
CONST lim = 1048576 '                 prime trial limit
DIM g AS STRING
CLS : tim# = TIMER: lin = 8

begin:
LOCATE 1, 1
PRINT "   +------------------------------------------------------+"
PRINT "   |  Input a,N to test if a belongs to phi(N) modulo N,  |"
PRINT "   |  input a = 1 to find the least primitive root of N,  |"
PRINT "   |  a = 0 to exit.                                      |"
PRINT "   |                                                      |"
PRINT "   +------------------------------------------------------+"
LOCATE 5, 7: INPUT ; "a"; g: fl = -1
IF INSTR(g, "'") THEN GOTO begin
IF VAL(g) < 1 THEN GOTO eind
IF VAL(g) = 1 THEN g = "2": fl = 0
INPUT " N"; Ni$
IF VAL(Ni$) < 2 THEN GOTO begin

LOCATE lin, 1
Readst a, g
Readst N, Ni$

IF fl THEN
   Gcd a, N, q
   Printn a, "a = ", "", 2
   IF NOT Isf(q, 1) THEN
      Printn N, "N = ", "", 2
      Prints "(a/N) = 0", 2: GOTO begin
   END IF
END IF

IF Isf(N, 2) THEN
   Prints "N = 2", 1
   IF fl THEN
      Prints "a is primitive Mod N", 2
   ELSE
      Prints "G = 1", 2
   END IF
   GOTO begin
ELSEIF Isf(N, 4) THEN
   Prints "N = 4", 1
   IF fl THEN
      IF Mp2(a, 4) = 3 THEN
         Prints "a is primitive Mod N", 2
      ELSE
         Prints "ord_N(a) = 1", 2
      END IF
   ELSE
      Prints "G = 3", 2
   END IF
   GOTO begin
END IF

s = Mp2(N, 4)
IF s = 2 THEN
   Rsft N, 1
   eveN = -1: Ni$ = "N = 2*"
ELSE
   eveN = 0: Ni$ = "N = "
END IF

sw = s > 0: ex = 0: g = ""
IF sw THEN
   Dup N, q
   DO
      IF IsPPrm(q, p, 5) THEN '       Miller-Rabin test
         EXIT DO
      ELSE
         s = 2: DO
            Letf b, CLNG(s) '         find quadratic residue b
            IF Kronec(b, q) = 1 THEN EXIT DO
            s = s + 1
         LOOP
         Modpwr b, q, q
         Inc b, -s: Gcd b, q, p
         IF Isf(p, 1) OR Cmp(p, q) = 0 THEN
            sw = 0: EXIT DO '         N is not a prime power
         END IF
         Swp p, q: ex = -1
      END IF
   LOOP
   '
   IF ex AND sw THEN
      Dup N, p: k = 0
      DO '                            divide out factor q
         Dup p, t0: Divd t0, q, t1
         IF NOT Isf(t0, 0) THEN EXIT DO
         Swp p, t1: k = k + 1
      LOOP
      IF NOT Isf(p, 1) THEN
         sw = 0
      ELSEIF k > 1 THEN
         Swp N, q '                   N:= N^ 1/k
         g = "^" + LTRIM$(STR$(k))
      END IF
   END IF
END IF

Printn N, Ni$, g, 2
IF sw = 0 THEN
   Prints "primitive root doesn't exist", 2
   GOTO begin
END IF

Dup N, R
Inc R, -1: Rsft R, 1 '                R:= (N - 1)/ 2
Dup R, b: k = 1 + Odd(b) '            N - 1 = 2^ k * odd(b)
g = "": IF k > 1 THEN g = "^" + LTRIM$(STR$(k))
Prints "    2" + g, 1

Triald q, b, lim '                    trial divide b
IF NOT Isf(b, 1) THEN
   sw = IsSqr(b, p)
   IF sw THEN Swp b, p
   IF IsPPrm(b, p, -5) THEN
      k = 1: IF sw THEN k = 2
      EnQ q, b, k
   ELSE
      Printn b, "fail: composite cofactor p-1 = F*R", "", 2
      Prints "", 1: GOTO begin
   END IF
END IF

Setl b, 0: t = 0
WHILE ExQ(p, k, q, t)
   g = "": IF k > 1 THEN g = "^" + LTRIM$(STR$(k))
   Printn p, "    ", g, 1
   Dup R, t0: Divd t0, p, t1 '        b:= set of R/ L, with
   EnQ b, t1, 1 '                     L odd primefactor of R
WEND
Swp q, b: Lsft R, 1

s = 2
DO '                                  find primitive root
   sw = Kronec(a, N) = -1 '           non-residue a
   IF eveN THEN sw = sw * -Mp2(a, 2) ' parity check
   IF sw THEN
      t = 0
      DO WHILE ExQ(p, 1, q, t) '      check all L[t]-Sylows
         Dup a, b
         Modpwr b, p, N
         IF Isf(b, 1) OR Ismin1(b, N) THEN
            sw = 0: EXIT DO
         END IF
      LOOP
      '
      IF ex AND sw THEN '             prime power check
         Dup a, b: Squ N, p
         Modpwr b, R, p
         IF Isf(b, 1) THEN sw = 0 '   ord_N^ 2(a) < phi(N^ 2)
      END IF
   END IF
   '
   IF fl THEN
      IF sw THEN
         g = "a is primitive Mod N"
      ELSE
         g = "ord_N(a) < phi(N)"
      END IF
      Prints g, 1: EXIT DO
   ELSE
      IF sw THEN '                    primitive root found
         Printn a, "G = ", "    ", 1
         EXIT DO
      ELSE '                          next a
         s = s + 1
         IF s = 4 OR s = 16 THEN '    skip perfect powers < 25
            s = s + 1
         ELSE
            IF s = 8 THEN s = 10
         END IF
         Letf a, CLNG(s): lin = CSRLIN
         PRINT "  a ="; s; : LOCATE lin, 1
      END IF
   END IF
LOOP
Prints "", 1
lin = CSRLIN
GOTO begin

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

