'******************************************************************************
'Subject: Multiplicative group of units modulo small N:
'         representing U(N) as a direct product of
'         canonical least-generator cyclic subgroups.
'Ref.   : D. Shanks, Solved and unsolved problems in number theory,
'         Chelsea, New York, 1978, pp. 83-98:
'        'Questions & answers concerning cycle graphs.'
'Author : Sjoerd.J.Schaper
'Date   : 07-03-2007
'Code   : QBasic, PDS, VBdos, FreeBasic (-lang qb)

'******************************************************************************
'This program is copyright (c) 2007 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'******************************************************************************
DEFINT A-Z
TYPE Pls
   f AS INTEGER
   e AS INTEGER
END TYPE
DECLARE FUNCTION gcd (a, b)
'highest common divisor
DECLARE SUB Pwr (a, k)
'let a:= a ^ k
DECLARE FUNCTION modpwr (a, k)
'return a ^ k (mod N)
DECLARE FUNCTION order (a)
'return ord_N(a)
DECLARE FUNCTION triald (a, sw)
'primefactors(a),
'set sw = 1 to print
DECLARE FUNCTION invars (a)
'torsion invariants(a)
DECLARE SUB Sieve (Ui, N)
'sieve residues prime to N
DECLARE SUB Balc (p)
'print balanced residues,
'initialize by printing N
DECLARE SUB Cycle (Hi, q, sw)
'H:= cycle generated by q,
'set sw = 1 to print
DECLARE SUB Smult (Ki, j, q)
'recursive set multiplication
DECLARE FUNCTION dprod (Gi, Hi)
'direct product G:= G x H,
'return 1 if G is enlarged
DECLARE SUB Remove (Gi, Hi)
'mask out subgroup H in G
DECLARE SUB Build (v())
'build generator, and cycle
DECLARE SUB Combi (v(), i)
'recursive index combinations
DECLARE FUNCTION popchk (Gi)
'return population count(Gi)

CONST U = 0, G = 1, H = 2      '       groups
CONST PB = 32, Pm1 = PB - 1    '       Dword size
DIM SHARED bi(Pm1) AS LONG     '       bitmasks
DIM SHARED bv(8192) AS LONG    '       4 bitvectors
DIM SHARED pt(3), ge(6), cy(6) '       generators, invariant orders
DIM SHARED pr(5) AS Pls        '       primelist
DIM SHARED cx, sx, fx, N       '       bounds, modulus
DIM iv(6)                      '       index vector
bi(0) = 1
FOR t = 1 TO Pm1 - 1 '                 store bitmasks
   bi(t) = bi(t - 1) * 2
NEXT t: bi(t) = &H80000000
CLS : tim# = TIMER

begin:
PRINT
DO: INPUT " N"; Ni$
LOOP WHILE INSTR(Ni$, "'")
IF VAL(Ni$) > 32767 THEN GOTO begin
N = VAL(Ni$): IF N < 2 THEN GOTO eind
Verb = RIGHT$(Ni$, 1) = "c" '          suffix 'c' to print cycles

Balc N: cx = invars(N): fi = 1 '       factorize N, find invariants
IF cx = 0 THEN PRINT " phi = 1": GOTO begin
cx = cx - 1
FOR c = 0 TO cx: fi = fi * cy(c): NEXT
PRINT " phi ="; fi '                   order(U(N))

sx = (N - 1) \ PB + 1
FOR r = 0 TO 3: pt(r) = r * sx: NEXT ' initialize row pointers
sx = sx - 1

bv(pt(G)) = 0: Sieve U, N '            prepare G and U(N)
fx = triald(cy(0), 0) - 1 '            factorize group exponent

FOR c = 0 TO cx '                      find generators
   FOR s = 0 TO sx
      IF bv(s) THEN '                  unpack U(N) elements
         r = s * PB
         FOR t = 0 TO Pm1
            IF bv(s) AND bi(t) THEN
               q = r + t
               IF order(q) = cy(c) THEN
                  Cycle H, q, 0 '      invariant order

                  IF dprod(G, H) THEN ' valid factor H
                     Remove U, G
                     ge(c) = q: GOTO jump
                  END IF
               END IF
            END IF
         NEXT t
      END IF
   NEXT s
jump:
NEXT c

IF popchk(G) = fi THEN '               check if G holds the full
   PRINT "<generator> order" '         unit group of order phi(N)
   FOR c = 0 TO cx
      PRINT " <"; ge(c); ">"; cy(c) '  minimal generating set
   NEXT c
   '
   IF Verb THEN
      PRINT "(indices) {cycle}  order"
      Combi iv(), 0 '                  covering cycles
      IF popchk(G) THEN
         PRINT " fail: units left"
      END IF
   END IF
ELSE
   PRINT " fail: order < phi"
END IF

GOTO begin

eind:
PRINT : PRINT "Timer:"; CSNG(TIMER - tim#); " s"
SYSTEM

SUB Balc (p)
STATIC d, r
DIM s AS STRING
   IF p < N THEN
      IF p > d THEN
         s = STR$(p - N)
      ELSE
         s = STR$(p)
      END IF
      PRINT SPACE$(r - LEN(s)); s; " ";
   ELSE
      d = N \ 2: r = LEN(STR$(d))
      PRINT " N ="; N
   END IF
END SUB

SUB Build (v())
   q = 1
   FOR c = 0 TO cx
      p = modpwr(ge(c), v(c))
      q = CLNG(p) * q MOD N
   NEXT c
   '
   p = pt(G) + q \ PB '                assuming G holds U(N)
   IF bv(p) AND bi(q AND Pm1) THEN
      PRINT " (";
      FOR c = 0 TO cx
         PRINT STR$(v(c));
      NEXT c: PRINT " )";
      Cycle H, q, 1
      Remove G, H '                    clear < q > candidates
   END IF
END SUB

SUB Combi (v(), c)
   FOR i = 1 TO cy(c)
      v(c) = i
      IF c < cx THEN
         Combi v(), c + 1
      ELSE
         Build v()
      END IF
   NEXT i
END SUB

SUB Cycle (r, q, sw)
   Hr = pt(r)
   FOR s = 0 TO sx: bv(Hr + s) = 0: NEXT
   '
   IF sw THEN k = 0: PRINT " {";
   p = 1: DO
      p = CLNG(q) * p MOD N
      s = Hr + p \ PB
      bv(s) = bv(s) OR bi(p AND Pm1) ' store cycle
      IF sw THEN k = k + 1: Balc p
   LOOP UNTIL p = 1
   IF sw THEN PRINT "} "; k
END SUB

FUNCTION dprod (i, j)
DIM W AS LONG
dprod = 0
   Gr = pt(i): Hr = pt(j): Kr = pt(j + 1)
   '
   IF bv(Gr) = 0 THEN
      SWAP pt(i), pt(j) '              1st largest subgroup H
      dprod = 1: EXIT FUNCTION
   END IF
   '
   bv(Hr) = bv(Hr) AND NOT bi(1) '     clear {1}
   FOR s = 0 TO sx
      W = bv(Gr + s) '                 non-trivial intersection:
      IF W AND bv(Hr + s) THEN EXIT FUNCTION
      bv(Kr + s) = W OR bv(Hr + s) '   join 1G with 1H
   NEXT s
   bv(Gr) = bv(Gr) AND NOT bi(1)
   '
   Smult j + 1, j, 0 '                 product G x H/1G, /1H
   SWAP pt(i), pt(j + 1)
dprod = 1
END FUNCTION

FUNCTION gcd (x, y)
   a = x: b = y
   WHILE b
      c = a MOD b
      a = b: b = c
   WEND
gcd = a
END FUNCTION

FUNCTION invars (a)
invars = 0
   fx = triald(a, 1) '                 factorize a
   IF fx = 0 THEN EXIT FUNCTION
   fx = fx - 1

   j = 0: t = -1
   IF pr(0).f = 2 THEN
      j = 1
      IF pr(0).e > 1 THEN '            multiplicity of 2
         e = pr(0).e: t = t + 1: cy(t) = 2
         IF e > 2 THEN t = t + 1: cy(t) = bi(e - 2)
      END IF
   END IF
   '
   FOR i = j TO fx '                   odd prime power decomposition
      q = 1: t = t + 1: cy(t) = pr(i).f
      IF pr(i).e > 1 THEN
         e = pr(i).e: q = cy(t)
         IF e > 2 THEN Pwr q, e - 1
         cy(t) = cy(t) * q
      END IF
      cy(t) = cy(t) - q '              Eulerphi(p^ e)
   NEXT i
   '
   FOR i = t TO 1 STEP -1 '            Smith canonical form
      FOR j = i - 1 TO 0 STEP -1
         IF cy(j) MOD cy(i) THEN
            q = gcd(cy(i), cy(j))
            cy(j) = (cy(j) \ q) * cy(i)
            cy(i) = q
         END IF
      NEXT j
   NEXT i
   '
   FOR i = t + 1 TO 6: cy(i) = 0: NEXT
invars = t + 1
END FUNCTION

FUNCTION modpwr (a, k)
DIM b AS LONG
   p = 1: b = a: e = k
   WHILE e
      IF e AND 1 THEN p = b * p MOD N
      b = b * b MOD N: e = e \ 2
   WEND
modpwr = p
END FUNCTION

FUNCTION order (a)
   e = cy(0) '                         Carmichael lambda
   FOR i = 0 TO fx
      FOR j = 1 TO pr(i).e
         k = e \ pr(i).f
         IF modpwr(a, k) = 1 THEN '    reduce exponent
            e = k
         ELSE
            EXIT FOR
         END IF
      NEXT j
   NEXT i
order = e
END FUNCTION

FUNCTION popchk (r)
DIM W AS LONG: Gr = pt(r)
   k = 0
   FOR s = 0 TO sx
      IF bv(Gr + s) THEN
         W = bv(Gr + s)
         FOR t = 0 TO Pm1
            k = k + (((W AND bi(t)) <> 0) AND 1)
         NEXT t
      END IF
   NEXT s
popchk = k
END FUNCTION

SUB Pwr (a, k)
   p = 1: e = k
   WHILE e
      IF e AND 1 THEN p = p * a
      a = a * a: e = e \ 2
   WEND: a = p
END SUB

SUB Remove (i, j)
   Gr = pt(i): Hr = pt(j)
   FOR s = 0 TO sx
      bv(Gr + s) = bv(Gr + s) AND NOT bv(Hr + s)
   NEXT s
END SUB

SUB Sieve (r, N)
DIM W AS LONG
   sw = pr(0).f = 2
   j = 0: W = &HFFFFFFFF
   IF sw THEN j = 1: W = &HAAAAAAAA
   FOR s = 0 TO sx: bv(s) = W: NEXT '  initialize bitvector U
   bv(r) = bv(r) AND NOT 3
   '
   FOR i = j TO fx '                   sieve units
      p = pr(i).f: dp = p
      IF sw THEN dp = p + p
      DO
         s = p \ PB
         bv(s) = bv(s) AND NOT bi(p AND Pm1)
         IF p > N - dp THEN EXIT DO
         p = p + dp
      LOOP
   NEXT i
END SUB

SUB Smult (i, j, q)
DIM W AS LONG: Gr = pt(j)
   FOR s = 0 TO sx
      IF bv(Gr + s) THEN '             unpack set elements
         W = bv(Gr + s): r = s * PB
         FOR t = 0 TO Pm1
            IF W AND bi(t) THEN
               IF j > 1 THEN
                  Smult i, j - 1, r + t
               ELSE '                  orbit of q under G
                  p = CLNG(r + t) * q MOD N
                  k = pt(i) + p \ PB
                  bv(k) = bv(k) OR bi(p AND Pm1)
               END IF
            END IF
         NEXT t
      END IF
   NEXT s
END SUB

FUNCTION triald (a, sw)
triald = 0
   IF a < 2 THEN EXIT FUNCTION
   f = 2: df = 1: fl = 0
   b = a: t = -1
   DO
      e = 0: DO
         q = b \ f '                   trial division
         IF b - q * f THEN EXIT DO
         b = q: e = e + 1
      LOOP
      IF e THEN
         t = t + 1 '                   store prime and exponent
         pr(t).f = f: pr(t).e = e
      ELSEIF f > q THEN '              f > sqrt(a)
         IF b > 1 THEN
            t = t + 1 '                store coprime
            pr(t).f = b: pr(t).e = 1
         END IF
         EXIT DO
      END IF
      IF fl THEN
         df = 6 - df: f = f + df '     sieve multiples of 2 and 3
      ELSE
         f = f + df: df = df * 2: fl = f = 5
      END IF
   LOOP
   '
   IF sw AND (t > 0 OR pr(0).e > 1) THEN
      PRINT "     ";
      FOR i = 0 TO t
         e = pr(i).e: s$ = ""
         IF e > 1 THEN s$ = "^" + LTRIM$(STR$(e))
         IF i < t THEN s$ = s$ + ","
         PRINT LTRIM$(STR$(pr(i).f)) + s$;
      NEXT i: PRINT
   END IF
triald = t + 1
END FUNCTION