'Subject: Roots of a quadratic congruence with prime modulus N,
'         the Shanks-Tonelli algorithm. For r^2 mod N = -1,
'         if N mod 4 = 1 also the Diophantine equation x^2 + y^2 = N
'         is solved with Cornacchia's algorithm.
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST N = 0, m = 1, r = 2, i1 = 3, i2 = 4
i = LargeInit(i2, "quadcngr")
DIM g AS STRING
CLS : tim# = TIMER: lin = 8

begin:
LOCATE 1, 1
PRINT "   +--------------------------------+"
PRINT "   |  input                         |"
PRINT "   |  a, N to solve r^2 mod N = a,  |"
PRINT "   |  0 to exit.                    |"
PRINT "   |                                |"
PRINT "   +--------------------------------+"
LOCATE 5, 7: INPUT ; "a"; ai$: fl = 0
IF INSTR(ai$, "'") THEN GOTO begin
IF VAL(ai$) = 0 THEN GOTO eind
INPUT ; " N"; Ni$
IF VAL(Ni$) < 2 THEN GOTO begin

LOCATE lin, 1
Readst i1, ai$
Readst N, Ni$: Sets N, 1
IF Gets(i1) THEN Moddiv i1, N
fl = Ismin1(i1, N) '                   a = N - 1: set flag to factorize N
IF fl THEN '                           in Gaussian primes (x + yi)(x - yi)
   Prints " a = N - 1", 1
ELSE
   Printn i1, " a = ", "", 2
END IF
Printn N, " N = ", "", 2

g = ""
IF Mp2(N, 2) = 0 THEN '               preliminary tests
   IF Isf(N, 2) THEN
      g = " r = 1  [1]"
      IF Mp2(i1, 2) = 0 THEN g = "(a/N)= 0"
   ELSE
      g = "N is even"
   END IF
ELSE
   t = Kronec(i1, N)
   IF t = 1 THEN
      t = IsPPrm(N, i2, 5)
      IF t = 0 THEN g = "N is composite"
   ELSE
      g = "(a/N)= -1"
      IF t = 0 THEN g = "(a/N)= 0"
   END IF
END IF
IF g <> "" THEN
   Prints g, 2: GOTO begin
END IF

Dup N, m
Inc m, -1: i = Odd(m) '                determine N - 1 = (2^ i) * odd(m)

Dup m, t3
Inc t3, 1: Rsft t3, 1 '                (m + 1)/ 2 in t3
Dup i1, r: Modpwr r, t3, N '           approximate root: r = a^(m+1)/2

IF i > 1 THEN '                        If i = 1 then N mod 4 = 3, root = r
   Euclid i1, N, t3 '                  a^ -1 (mod N) in i1
   '
   b = 2 '                             Shanks-Tonelli algorithm Ressol:
   DO
      Letf i2, CLNG(b)
      IF Kronec(i2, N) = -1 THEN '     find non-residue b
         Modpwr i2, m, N: EXIT DO '    generator z = b^ m (mod N) in i2
      END IF
      b = b + 1
   LOOP
   '
   t = 0: q = i
   DO
      Squ r, t0: Divd t0, N, t1
      Modmlt t0, i1, N '               t0:= r^ 2 * a^ -1 mod N
      p = q: q = 1
      DO '                             find exponent
         IF Isf(t0, 1) THEN
            IF q = 1 THEN GOTO jump '  root found
            EXIT DO '                  ord_N(t0) = 2^ q
         END IF
         Modsqu t0, N: q = q + 1
      LOOP
      FOR i = q TO p - 1 '             reduce exponent
         Modsqu i2, N
      NEXT i
      Modmlt r, i2, N: t = t + 1 '     r:= r * z^ (2^ p-q) mod N
   LOOP
jump:
   Prints " b =" + STR$(b) + "  t =" + STR$(t), 1
END IF
Modbal r, N: Sets r, 1
Printn r, " r = ", "", 2

IF fl THEN
   Isqrt N, i1 '                       Cornacchia's algorithm calculates
   Dup N, i2 '                         the quadratic representation of N.
   t = 0: DO
      t = t + 1 '                      Euclidean division (N, root)
      Divd i2, r, m '                  quotients m form a symmetric,
      IF Isf(i2, 0) THEN '             finite continued fraction.
         IF t AND 1 THEN '             make even cf length
            Inc m, -1
            Printn m, "", ",", 0
            t = t + 1: Letf m, 1
         END IF
         Printn m, "", "  [" + LTRIM$(STR$(t)) + "]", 1
         EXIT DO
      END IF
      IF fl THEN
         IF Cmp(i2, i1) < 1 THEN '     remainder <= Isqrt(N)
            fl = 0: Dup i2, i1 '       store remainder = x, skip search
         END IF
      END IF
      Swp i2, r
      Printn m, "", ",", 0
   LOOP
   Printn i1, " x = ", "", 2
   Squ i1, i2: Subt N, i2: Isqrt N, i2
   Printn i2, " y = ", "", 2 '         y = Isqrt(N - x^ 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

