'Subject: Prime factorization and arithmetic functions.
'         Trial division: quick check for factors up to 2^ 21.
'         Calculates Euler totient and Carmichael lambda functions,
'         also the divisor functions æ(a), w(a), ê(a), ë(a) and å(a).
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DECLARE SUB Afuncs (fr, k, sw)
CONST a = 0, dr = 1, q = 2, p = 3, l = 4, d = 5, s = 6
r = LargeInit(s, "trialdiv")
CONST lim = 2097143 '                 search limit
DIM SHARED sq
CLS : tim# = TIMER: lin = 7

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

LOCATE lin, 1
Readst a, ai$
Printn a, "    a = ", "", 2
Afuncs a, 0, 0
pr = 0: sq = 0
IF Getl(a) > 2 THEN '                 preliminary tests
   sq = IsSqr(a, q)
   IF sq THEN Swp a, q '              a is square
   pr = IsPPrm(a, q, 5) '             a or sqrt(a) is pprime
END IF

IF pr THEN
   Afuncs a, 1, 1
ELSE
   Triald q, a, lim: t = 0
   WHILE ExQ(dr, k, q, t)
      Afuncs dr, k, 1
   WEND
   '
   IF NOT Isf(a, 1) THEN
      IF IsPPrm(a, q, -5) THEN
         Afuncs a, 1, 1 '             prime cofactor
      ELSE
         Prints " failure ", 2 '      quit job
         GOTO begin
      END IF
   END IF
END IF
Readst a, ai$
Afuncs a, 0, 2 '                      print results

Prints "", 1
lin = CSRLIN
GOTO begin

eind:
LOCATE lin, 1
g$ = STR$(CSNG(TIMER - tim#))
Prints "Timer:" + g$ + " s", 1
CALL Term
SYSTEM

SUB Afuncs (fr, k, sw)
STATIC bo, o

SELECT CASE sw
CASE 0 '                             initialize
   Letf p, 1: Letf l, 1: Letf s, 1 '  totient, lambda & divisor functions
   Letf d, 1: o = 0: bo = 0
   '
CASE 1 '                              main
   g$ = "": IF sq THEN k = k * 2
   IF k > 1 THEN g$ = "^ " + LTRIM$(STR$(k))
   Printn fr, SPACE$(8), g$, 1
   '
   o = o + 1: bo = bo + k '           omega functions
   Letf t0, k + 1: Lmul d, t0 '       s_0(a)
   Letf t3, 1: Letf t1, 1
   FOR i = 1 TO k
      Mult t3, fr, t0
      Add t1, t0: Swp t3, t0
   NEXT i
   Lmul s, t1 '                       s_1(a)
   Subt t3, t0: Lmul p, t3 '          phi(a)
   IF Isf(fr, 2) AND k > 2 THEN Rsft t3, 1
   Lcm l, t3, fr: Swp l, fr '         lambda(a)
   '
CASE ELSE
   Printn p, "totient  ", "", 2 '     size of the unit group (mod a)
   Printn l, " lambda  ", "", 2 '     largest cyclic subgroup (mod a)
   m = 1 - 2 * (o AND 1) '            Moebius function
   IF bo > o THEN m = 0
   Prints "moebius " + STR$(m), 0 '   parity,
   Prints "   omega " + STR$(o), 0 '   distinct primefactors
   Prints "   bigOo " + STR$(bo), 1 '   total primefactors
   Printn d, "numdivr  ", "", 1 '     number of divisors
   Printn s, "sumdivr  ", "", 2 '      sum of all divisors
   Lsft fr, 1: Subt s, fr
   Printn s, "aliquot  ", "", 2 '     abundance s - 2a
END SELECT
END SUB

