A unitary divisor of a number n is a number d such that d|n and gcd(d, n/d)=1. For example, 3 is a unitary divisor of 12 because gcd(3, 12/3)=gcd(3, 4)=1.

A unitary perfect number is a number that is the sum of its unitary divisors less than itself. For example, 60 is a unitary perfect number because its unitary divisors less than itself are 1, 3, 4, 5, 12, 15, and 20, and 1+3+4+5+12+15+20=60. Find all unitary perfect numbers less than 1000000.

The key to figuring this out (more satisfying than looking it up) is recognizing the need for the unitary divisor to have, for each prime factor in common with the original number (as each prime factor of the divisor must be), to be raised to the same power as in the original number or be raised to the zero power, that is, not to be a factor at all.

The program:

DECLARE SUB factor (num#)
DECLARE SUB addOn (which#)
DEFDBL A-Z
CLEAR , , 25000

OPEN "unitary perfect.txt" FOR OUTPUT AS #2

DIM SHARED fct(10, 2), nprmfact, unitdiv, sum, whole, ufct(1024), uct

FOR n = 2 TO 1000000
    unitdiv = 1: sum = 0: whole = n: uct = 0
    factor n
    addOn 1
    IF sum = n THEN
        PRINT #2, n
        FOR i = 1 TO nprmfact: PRINT #2, fct(i, 1); fct(i, 2): NEXT
        FOR i = 1 TO uct: PRINT #2, ufct(i);: NEXT
        PRINT #2,: PRINT #2,
    END IF
NEXT n

SUB addOn (which)
FOR incl = 0 TO 1
    IF incl THEN unitdiv = unitdiv * INT(fct(which, 1) ^ fct(which, 2) + .5)
    IF which = nprmfact THEN
        IF unitdiv < whole THEN
            uct = uct + 1
            ufct(uct) = unitdiv
            sum = sum + unitdiv
        END IF
    ELSE
        addOn which + 1
    END IF
    IF incl THEN unitdiv = unitdiv / INT(fct(which, 1) ^ fct(which, 2) + .5)
NEXT
END SUB

SUB factor (num)
nprmfact = 0: n = ABS(num): IF n > 0 THEN limit = SQR(n): ELSE limit = 0
IF limit <> INT(limit) THEN limit = INT(limit + 1)
dv = 2: GOSUB DivideIt
dv = 3: GOSUB DivideIt
dv = 5: GOSUB DivideIt
dv = 7
DO UNTIL dv > limit
    GOSUB DivideIt: dv = dv + 4 '11
    GOSUB DivideIt: dv = dv + 2 '13
    GOSUB DivideIt: dv = dv + 4 '17
    GOSUB DivideIt: dv = dv + 2 '19
    GOSUB DivideIt: dv = dv + 4 '23
    GOSUB DivideIt: dv = dv + 6 '29
    GOSUB DivideIt: dv = dv + 2 '31
    GOSUB DivideIt: dv = dv + 6 '37
    IF INKEY$ = CHR$(27) THEN s$ = CHR$(27): EXIT SUB
LOOP
IF n > 1 THEN nprmfact = nprmfact + 1: fct(nprmfact, 1) = n: fct(nprmfact, 2) = 1
EXIT SUB

DivideIt:
count = 0
DO
    q = INT(n / dv)
    IF q * dv = n AND n > 0 THEN
        n = q: count = count + 1: IF n > 0 THEN limit = SQR(n): ELSE limit = 0
        IF limit <> INT(limit) THEN limit = INT(limit + 1)
    ELSE
        EXIT DO
    END IF
LOOP
IF count > 0 THEN
    nprmfact = nprmfact + 1
    fct(nprmfact, 1) = dv
    fct(nprmfact, 2) = count
END IF
RETURN
END SUB

finds

 6
 2  1
 3  1
 1  3  2

 60
 2  2
 3  1
 5  1
 1  5  3  15  4  20  12

 90
 2  1
 3  2
 5  1
 1  5  9  45  2  10  18

 87360
 2  6
 3  1
 5  1
 7  1
 13  1
 1  13  7  91  5  65  35  455  3  39  21  273  15  195  105  1365  64  832  448  5824  320  4160  2240  29120  192  2496  1344  17472  960  12480  6720

where the first line in each group is the unitary perfect number itself; the next group of lines consist of each prime factor followed by the power to which it is raised in that number; and the last line (or group of lines in the last instance) constitute the number's unitary divisors that are less than the number itself. The last case shown, 87,360, has 2^5 - 1 = 31 unitary divisors shown as each of the prime factors 2, 3, 5, 7 and 13, can have its respective power either included or excluded; when all are included, it's the number itself and is not shown.

 

Posted by Charlie on 2013-12-16 16:26:07