'Subject: Large integer exponentiation [modulo N], 10,916 digits max.
'         The Miller-Rabin strong pseudoprime-test is performed for
'         exponent "N-1". Base a is inverted modulo N for exponent k < 0.
'         Input N = 0 to disable modulo reduction.
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST a = 0, k = 1, N = 2, d = 3
r = LargeInit(d, "powrmodn")
CLS : tim# = TIMER: lin = 8

begin:
LOCATE 1, 1
PRINT "   +-------------------------------------------------------+"
PRINT "   |  Input a, k, (N) to calculate a ^ k (modulo N)        |"
PRINT "   |  Input k = 'N-1' to perform strong pseudoprime-test,  |"
PRINT "   |  k =-1 to invert a mod N,  a = 0 to exit.             |"
PRINT "   |                                                       |"
PRINT "   +-------------------------------------------------------+"
LOCATE 5, 7: INPUT "a"; ai$
IF INSTR(ai$, "'") THEN '             comment string
   IF INSTR(LCASE$(ai$), "verb") THEN ' verbose mode
      Verb = -1
   ELSEIF LEN(ai$) > 1 THEN
      Verb = 0: Prints ai$, 1
   END IF
   GOTO begin
END IF
IF VAL(ai$) = 0 THEN GOTO eind
Readst a, ai$ '                       base a
'
LOCATE 5, 7: PRINT SPACE$(53) + "|"
LOCATE 5, 7: INPUT "k"; g$
fl = -SGN(INSTR(UCASE$(g$), "N")) '   set flag to apply primality testing
IF VAL(g$) < 0 THEN g$ = "-1" '       compute a^ -1 mod N
'
LOCATE 5, 7: PRINT SPACE$(53) + "|"
LOCATE 5, 7: INPUT "N"; Ni$
Readst N, Ni$: Sets N, 1 '            modulus abs(N)
sw = NOT Isf(N, 0) '                  set switch to apply reduction Mod N
IF fl AND sw THEN
   Dup N, k
ELSE
   Readst k, g$ '                     exponent k
END IF

LOCATE lin, 1
Printn a, "    a = ", "", 2
IF fl THEN
   Inc k, -1
   Prints "    k = N - 1", 1
ELSE
   IF sw THEN fl = Ismin1(k, N) '     free test
   Printn k, "    k = ", "", 2
END IF
IF sw THEN Printn N, "    N = ", "", 2

IF Isf(k, -1) THEN
   Euclid a, N, d '                   multiplicative inverse
   IF Isf(d, 1) THEN '                gcd = 1
      Printn a, "a^-1 ð ", " mod N", 1
   ELSE
      Prints "impossible inverse modulo N,", 1
      Printn d, "(a,N) = ", "", 1
   END IF
   '
ELSE
   IF sw THEN
      Gcd a, N, d: g$ = ""
      IF NOT Isf(d, 1) THEN g$ = "(a/N) = 0"
      IF fl AND Mp2(N, 2) = 0 THEN g$ = "  N is even"
      IF g$ <> "" THEN
         IF fl THEN g$ = g$ + ", testing cancelled"
         fl = 0: Prints g$, 1
      END IF
   END IF

   i = Odd(k) '                       determine k = (2^ i) * odd(m)
   Modpwr a, k, N '                   Miller-Rabin test: if N is composite
   Nc = NOT Isf(a, 1) '               then a^ m Mod N <> 1...
   FOR j = 1 TO i
      IF fl THEN
         r = Ismin1(a, N) '           and (a^ m)^ (2^ j) Mod N <> -1
         Nc = Nc AND NOT r '          for all j in [0, i - 1].
         nr = (j = i) AND r '         Non-residue if a^ (N-1)/2 Mod N = -1
      END IF

      IF Verb THEN
         Modbal a, N
         Printn a, "(a^m)^ 2^" + LTRIM$(STR$(j - 1)) + " = ", " mod N", 1
      END IF
      Squ a, t0: Swp a, t0
      IF sw THEN Divd a, N, t0
   NEXT j

   IF sw THEN
      Modbal a, N '                   output power
      Printn a, " a^ k ð ", " mod N", 1
   ELSE
      Printn a, " a^ k = ", "", 2
   END IF

   IF fl THEN '                       evaluate primality test
      IF Nc THEN
         g$ = "   N is composite"
      ELSE
         g$ = "   N is a strong pseudoprime to a, "
         h$ = "(a/N) = 1": IF nr THEN h$ = "(a/N) = -1"
         g$ = g$ + h$ '               Legendre's symbol
      END IF
      Prints g$, 1 '                  output test
   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

