'Subject: Shanks's square form factorization with queue and multipliers.
'         Non-trivial ambiguous forms of discriminant 4*kN yield factors of N,
'         fast, but only useful for small numbers with up to 20-25 digits.
'Author : Sjoerd.J.Schaper
'Date   : 31-01-2004
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DEFINT C
DECLARE SUB Rho (u, v, w, sw)
'reduce indefinite form, set sw = 0 if u is null
DECLARE FUNCTION proper (v, s)
'return -1 if a proper square form is found
DECLARE SUB Qform (g AS STRING, t AS LONG, u, v, w)
'print binary quadratic form(t) (u, 2v, w)
DECLARE SUB Queue (Q)
'print {Q}

CONST lim = 11, ix = 5000000 '         small prime- and search limits
CONST N = 0, a = 1, b = 2, c = 3
CONST r = 4, x = 5, y = 6, z = 7 '     largeint pointers
CONST kN = 8, rN = 9, cb = 10, rb = 11
CONST k = 12, Qb = 13, Qr = 14, p = 15
t = LargeInit(p, "squfofac")

DIM SHARED lin
DIM i AS LONG, j AS LONG
DIM mp(3): mp(0) = 1
mp(1) = 3: mp(2) = 5: mp(3) = 7 '      custom multipliers
CLS : tim# = TIMER: lin = 7

begin:
LOCATE 1, 1
PRINT "   +---------------------+"
PRINT "   |  input N to split,  |"
PRINT "   |  0 to exit.         |"
PRINT "   |                     |"
PRINT "   +---------------------+"
LOCATE 4, 7: INPUT ; "N"; Ni$
IF INSTR(Ni$, "'") THEN GOTO begin
IF VAL(Ni$) < 2 THEN GOTO eind

LOCATE lin, 1
Readst N, Ni$
Printn N, " N = ", "", 2

Triald Qb, N, lim: t = 0 '             remove small primefactors
IF Getl(Qb) THEN Prints " f =", 0
WHILE ExQ(p, e, Qb, t)
   g$ = "": IF e > 1 THEN g$ = "^" + LTRIM$(STR$(e))
   Printn p, " ", g$ + ",", 0
WEND
IF t THEN
   IF Isf(N, 1) THEN
      Prints "", 2
      lin = CSRLIN: GOTO begin
   ELSE
      Printn N, ".. N/f = ", "", 2
   END IF
END IF

IF IsPPrm(N, x, 5) THEN '              preliminary tests
   g$ = " N": IF t THEN g$ = " N/f"
   Prints g$ + " is probably prime", 2
   lin = CSRLIN: GOTO begin
ELSEIF IsSqr(N, x) THEN
   Printn x, "   = ", "^2", 2
   Prints "", 1
   lin = CSRLIN: GOTO begin
END IF

fl = 0
FOR tl = 0 TO 3 '                      main
   IF mp(tl) <= 1 THEN
      Dup N, kN: Letf k, 1
   ELSE
      Letf k, CLNG(mp(tl))
      Mult k, N, kN '                  kN = multiplier * N
      Printn k, " k = ", " ", 0
   END IF
   Isqrt kN, rN
   Printn rN, " r = ", "", 1 '         rN = [sqrt(kN)]

   Dup rN, b: Letf c, 1
   Rho a, b, c, 0
   Qform "P", 1, a, b, c '             principal form

   Setl Qb, 0: Setl Qr, 0
   Dup rN, r: Lsft r, 1
   Isqrt r, rb: Lsft k, 1 '            Q-bound 2k[sqrt(2rN)]
   Mult rb, k, cb

   FOR i = 2 TO ix '                   search principal cycle
      IF Cmp(c, cb) < 1 THEN
         Gcd c, k, p: Dup c, r
         IF NOT Isf(p, 1) THEN Ldiv r, p
         IF Cmp(r, rb) < 1 THEN '      enqueue b and c/gcd(c,2k)
            EnQ Qb, b, 0: EnQ Qr, r, 0
         END IF
      END IF

      Rho a, b, c, -1
      IF (i AND 1) = 0 THEN
         IF (i AND 32767) = 0 THEN
            LOCATE lin, 5: PRINT i;
         END IF

         IF IsSqr(c, r) THEN '         square form found
            IF Isf(c, 1) THEN EXIT FOR
            '
            IF proper(b, r) THEN
               Qform "P", i, a, b, c

               Gcd b, r, p: j = 1
               IF Isf(p, 1) THEN
                  Dup b, y: Swp r, z
                  Chs y: Rho x, y, z, 0
                  Qform "A", j, x, y, z ' inverse square-root
                  '
                  DO
                     j = j + 1 '       search ambiguous cycle
                     IF (j AND 32767) = 0 THEN
                        LOCATE lin, 5: PRINT j;
                     END IF

                     Dup y, r
                     Rho x, y, z, -1
                     IF Cmp(y, r) = 0 THEN '   symmetry point
                        '
                        Gcd N, x, p
                        IF Isf(p, 1) THEN EXIT DO
                        Qform "A", j, x, y, z
                        Swp x, p: fl = -1: GOTO jump
                     END IF
                  LOOP
 
               ELSE '                  gcd(b,sqrt(c)) > 1
                  Squ p, r: Gcd N, r, x
                  fl = -1: GOTO jump
               END IF
            END IF
         END IF
      END IF
   NEXT i: PRINT
NEXT tl

jump:
IF fl THEN
   g$ = " f = "
   IF i THEN
      IF Getl(Qr) THEN Queue Qr
      IF j = 1 THEN g$ = "gcd^2 = "
   END IF
   Printn x, g$, " ", 0
   Divd N, x, y
   Printn y, " N/f = ", "", 2
ELSE
   Prints " fail", 1
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

FUNCTION proper (v, s)
sw = -1: t = 0: Q = Qr: i = 0
DO WHILE ExQ(x, 0, Q, i)
   IF Cmp(s, x) = 0 THEN '             root s in queue
      ib = 0
      FOR j = 0 TO t
         u = ExQ(y, 0, Qb, ib)
      NEXT j: Subt y, v
      Divd y, s, t0
      '
      IF Isf(y, 0) THEN '     v = y mod s: conjugate factors
         FOR j = 0 TO 1
            Setl Q, Getl(Q) + 1
            Rsft Q, -i: i = ib
            Setl Q, Getl(Q) - 1: Q = Qb
         NEXT j
         sw = 0: EXIT DO '             update queue and skip
      END IF
   END IF
   t = t + 1
LOOP
proper = sw
END FUNCTION

SUB Qform (g AS STRING, t AS LONG, u, v, w)
IF t > 32767 THEN PRINT
IF t AND 1 THEN
   Chs w
ELSE
   Chs u
END IF
   g = g + LTRIM$(STR$(t)) + " = ("
   Printn u, g, ",", 0
   Lsft v, 1
   Printn v, "", ",", 0
   Rsft v, 1
   Printn w, "", ")", 1
   lin = CSRLIN
IF t AND 1 THEN
   Chs w
ELSE
   Chs u
END IF
END SUB

SUB Queue (Q)
   Prints "{Q} [", 0
   g$ = ",": i = 0
   WHILE ExQ(p, 0, Q, i)
      IF i = Getl(Q) THEN g$ = ""
      Printn p, "", g$, 0
   WEND
   Prints "]", 1
END SUB

SUB Rho (u, v, w, sw)
   Dup v, t0: Swp u, w '               partial quotient
   Add v, rN: Divd v, u, p '           p:= (B + sqrt(kN)) / C
   Subt v, rN: Chs v '                 B':= p * C - B
   IF sw THEN
      Subt t0, v
      Mult t0, p, t1: Add w, t1 '      fast residue w/previous C
   ELSE
      Squ v, t0: Subt t0, kN '         quadratic residue
      Chs t0: Divd t0, u, w '          C':= (kN - B'²) / C
   END IF
END SUB

