'******************************************************************************
'Subject: Prime factorization for small Gaussian integers
' with Norm(a + bi) < 2^ 31.
' Calculates Euler totient and Carmichael lambda functions,
' also the divisor functions d(c), s(c), o(c), O(c) and mu(c).
'Author : Sjoerd.J.Schaper
'Date : 08-14-2006
'Code : FreeBasic and all QBasic's
'You'll find a version for big integer-BASIC
here
'******************************************************************************
'This program is copyright (c) 2006 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.
'******************************************************************************
REM $LANG: "qb"
DEFLNG A-Z
TYPE complex
r AS DOUBLE
i AS DOUBLE
END TYPE
DECLARE SUB triald (N)
'factorize norm(N) by trial division
DECLARE SUB afuncs (p, sw)
'arithmetical functions
DECLARE SUB quadrep (a AS complex, p)
'quadratic representation of p
DECLARE FUNCTION modpow# (a, b, c)
'modular exponentiation
DECLARE FUNCTION modd# (a AS DOUBLE, b)
'floating point a modulo b
DECLARE FUNCTION lcm (a, b)
'least common multiple
DECLARE SUB cadd (a AS complex, b AS complex)
DECLARE SUB csub (a AS complex, b AS complex)
DECLARE SUB cmul (a AS complex, b AS complex)
DECLARE SUB cgcd (a AS complex, b AS complex)
'complex arithmetic, results are returned in a
DECLARE SUB cdivmod (q AS complex, r AS complex, a AS complex, b AS complex)
'quotient(a/b) is returned in q, remainder in r
DECLARE SUB unit (a AS complex)
DECLARE FUNCTION norm# (a AS complex)
DECLARE FUNCTION cprn$ (a AS complex)
'pretty print complex number
'
CONST Eps = 1E-08
CONST mxd# = 9007199254740991# ' full 53 bit mantisse
DIM SHARED errsw AS INTEGER ' overflow in modd()
DIM SHARED Nc AS complex ' the number
CLS : LOCATE 2, 1: tim# = TIMER
'
DO
INPUT " re(c) ", Nc.r
INPUT " im(c) ", Nc.i
IF norm(Nc) > &H7FFFFFFF THEN
PRINT " norm(N) too big"
ELSE
N = norm(Nc)
IF N <= 1 THEN EXIT DO
PRINT " C = "; cprn(Nc)
PRINT " norm "; N
errsw = 0: triald N
END IF
PRINT
LOOP
PRINT " time"; TIMER - tim#
SYSTEM
SUB afuncs (N, sw) STATIC
DIM a AS complex, dr AS complex, rm AS complex
DIM s AS complex, qt AS complex
'
SELECT CASE sw
CASE 0 ' initialize
p = 1: l = 1: unit s
d = 1: o = 0: bo = 0
a = Nc
'
CASE 1 ' main
IF (N AND 3) = 3 THEN
dr.r = N: dr.i = 0
ph = N * N: tx = 0
ELSE
quadrep dr, N ' split N = dr^ 2 + di^ 2
IF errsw THEN EXIT SUB
ph = N: tx = ABS(N > 2) ' 0 / 1
END IF
'
FOR t = 0 TO tx
k = 0
DO
cdivmod qt, rm, a, dr ' trial divide a
IF norm(rm) > Eps THEN EXIT DO
SWAP a, qt: k = k + 1 ' running quotient a:= a/dr
LOOP
IF k > 0 THEN ' factor found
g$ = ""
IF k > 1 THEN g$ = " ^" + LTRIM$(STR$(k))
PRINT SPACE$(8); cprn(dr); g$
o = o + 1: bo = bo + k ' omega functions
'
t0 = 1
FOR i = 1 TO k
t1 = t0
t0 = t0 * ph
NEXT i
t0 = t0 - t1: p = p * t0 ' phi(a) p^ k - p^(k-1)
IF ph = 2 AND k > 2 THEN t0 = t0 \ 2
l = lcm(l, t0) ' lambda(Nc)
'
d = d * (k + 1) ' s_0(Nc)
unit rm: unit qt
FOR i = 1 TO k
cmul rm, dr
cadd qt, rm ' 1 + p +... + p^ k
NEXT i
cmul s, qt ' s_1(Nc)
END IF
SWAP dr.r, dr.i ' conjugate(dr)
NEXT t
'
CASE ELSE
PRINT " totient "; p ' size of the unit group (mod Nc)
PRINT " lambda "; l ' largest cyclic subgroup (mod Nc)
m = 1 - 2 * (o AND 1)
IF o < bo THEN m = 0 ' Moebius function
PRINT " moebius "; m; " omega "; o; " bigOo "; bo
PRINT " numdivr "; d ' number of divisors
PRINT " sumdivr "; cprn(s) ' sum of all divisors
csub s, Nc
PRINT " aliquot "; cprn(s) ' sum of proper divisors
PRINT " norm "; norm(s)
END SELECT
END SUB
SUB cadd (a AS complex, b AS complex)
a.r = a.r + b.r
a.i = a.i + b.i
END SUB
SUB cdivmod (q AS complex, r AS complex, a AS complex, b AS complex)
DIM d AS DOUBLE: d = norm(b)
q = a: b.i = -b.i
cmul q, b: b.i = -b.i
q.r = INT(q.r / d + .5)
q.i = INT(q.i / d + .5)
r = q: cmul r, b
r.r = a.r - r.r
r.i = a.i - r.i
END SUB
SUB cgcd (a AS complex, b AS complex)
DIM q AS complex, r AS complex
WHILE norm(b) > Eps
cdivmod q, r, a, b
a = b: b = r
WEND
END SUB
SUB cmul (a AS complex, b AS complex)
DIM tmp AS DOUBLE
tmp = a.r * b.r - b.i * a.i
a.i = a.i * b.r + b.i * a.r
a.r = tmp
END SUB
FUNCTION cprn$ (a AS complex)
DIM g AS STRING, s AS STRING
g = STR$(a.r)
IF ABS(a.i) > Eps THEN
s = " + ": IF a.i < 0 THEN s = " - "
IF ABS(a.i) = 1 THEN
g = g + s + "i"
ELSE
g = g + s + LTRIM$(STR$(ABS(a.i))) + "i"
END IF
END IF
cprn = g
END FUNCTION
SUB csub (a AS complex, b AS complex)
a.r = a.r - b.r
a.i = a.i - b.i
END SUB
FUNCTION lcm (a, b)
u = a: v = b
WHILE v > 0
t = v
v = u MOD v
u = t
WEND
lcm = a * (b \ u)
END FUNCTION
FUNCTION modd# (a AS DOUBLE, b)
errsw = a > mxd ' floating point gap > .5
modd = a - INT(a / b) * b
END FUNCTION
FUNCTION modpow# (a, b, c)
DIM p AS DOUBLE, s AS DOUBLE
p = 1: s = a: k = b
WHILE k ' R=>L binary powering
IF k AND 1 THEN
p = modd(p * s, c)
END IF
s = modd(s * s, c)
k = k \ 2
WEND
modpow = p
END FUNCTION
FUNCTION norm# (a AS complex)
norm = a.r * a.r + a.i * a.i
END FUNCTION
SUB quadrep (a AS complex, N)
DIM r AS complex
IF N = 2 THEN
a.r = 1: a.i = 1
ELSE
N1 = N - 1: N2 = N1 \ 2 ' solve r^ 2 = -1 (mod N)
FOR b = 2 TO 999 ' find non-residue b
IF modpow(b, N2, N) = N1 THEN EXIT FOR
IF errsw THEN EXIT SUB
NEXT b
r.r = modpow(b, N2 \ 2, N) ' root = b^ (N - 1)/ 4 (mod N)
a.r = N: a.i = 0 ' N divides r^ 2 + 1 = (r + i)(r - i)
r.i = 1: cgcd a, r ' a = gcd(N, r + i)
a.r = ABS(a.r): a.i = ABS(a.i)
IF a.r < a.i THEN SWAP a.r, a.i
END IF
END SUB
SUB triald (N)
afuncs N, 0 ' initialize
'
a = N: sw = 0
dr = 2: df = 1
DO
k = -1
DO
t1 = a \ dr ' trial divide a
t0 = a - t1 * dr ' remainder
IF t0 THEN EXIT DO
SWAP a, t1: k = 0 ' running quotient a:= a/dr
LOOP
IF k THEN
IF dr > t1 THEN EXIT DO ' divisor > sqrt(a)
ELSE
afuncs dr, 1 ' factor found
IF errsw THEN EXIT DO
END IF
IF sw THEN
DO
df = 6 - df: dr = dr + df ' sieve multiples of 2, 3 and 5
LOOP UNTIL dr MOD 5
ELSE
dr = dr + df: df = df * 2
sw = dr = 5
END IF
LOOP
'
IF a > 1 THEN afuncs a, 1 ' nontrivial factor a
IF errsw THEN
PRINT " overflow"
ELSE
afuncs a, 2 ' print results
END IF
END SUB
SUB unit (a AS complex)
a.r = 1: a.i = 0
END SUB