Practical numbers
A Practical number P has some selection of its proper divisors, (other than itself), that can be selected to sum to every integer less than itself.
Compute all the proper divisors/factors of an input number X, then, using all selections from the factors compute all possible sums of factors and see if all numbers from 1 to X-1 can be created from it.
- Task
Write a function that given X returns a boolean value of whether X is a Practical number, (using the above method).
- Show how many Practical numbers there are in the range 1..333, inclusive.
- Show that the Practical numbers in the above range start and end in:
- 1, 2, 4, 6, 8, 12, 16, 18, 20, 24 ... 288, 294, 300, 304, 306, 308, 312, 320, 324, 330
- Stretch Goal
- Show if 666 is a Practical number
Phix
atom t0 = time() function sumset(sequence res, set, integer tot=0, done=0) if tot<length(res) then if tot!=0 and res[tot]=0 then res[tot] = 1 res[$] += 1 end if done += 1 if done<=length(set) and res[$]<length(res) then res = sumset(res,set,tot+set[done],done) res = sumset(res,set,tot,done) end if end if return res end function function is_practical(integer n) return not find(0,sumset(repeat(0,n),factors(n,-1))) end function sequence res = apply(true,sprintf,{{"%3d"},filter(tagset(333),is_practical)}) printf(1,"Found %d practical numbers:\n%s\n\n",{length(res),join(shorten(res,"",10),", ")}) procedure stretch(integer n) printf(1,"is_practical(%d):%t (%s)\n",{n,is_practical(n),elapsed(time()-t0)}) end procedure papply({666,6666,66666,672,720},stretch)
- Output:
Found 76 practical numbers: 2, 4, 6, 8, 12, 16, 18, 20, 24, 28, ..., 288, 294, 300, 304, 306, 308, 312, 320, 324, 330 is_practical(666):true (0.3s) is_practical(6666):true (0.3s) is_practical(66666):false (0.3s) is_practical(672):true (1.5s) is_practical(720):true (1 minute and 7s)
Python
Python: Straight forward implementation
<lang python>from itertools import chain, cycle, accumulate, combinations from typing import List, Tuple
- %% Factors
def factors5(n: int) -> List[int]:
"""Factors of n, (but not n)."""
def prime_powers(n):
# c goes through 2, 3, 5, then the infinite (6n+1, 6n+5) series
for c in accumulate(chain([2, 1, 2], cycle([2,4]))):
if c*c > n: break
if n%c: continue
d,p = (), c
while not n%c:
n,p,d = n//c, p*c, d + (p,)
yield(d)
if n > 1: yield((n,))
r = [1]
for e in prime_powers(n):
r += [a*b for a in r for b in e]
return r[:-1]
- %% Powerset
def powerset(s: List[int]) -> List[Tuple[int, ...]]:
"""powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3) .""" return chain.from_iterable(combinations(s, r) for r in range(1, len(s)+1))
- %% Practical number
def is_practical(x: int) -> bool:
""" Is x a practical number.
I.e. Can some selection of the proper divisors of x, (other than x), sum
to i for all i in the range 1..x-1 inclusive.
"""
if x == 1:
return True
if x %2:
return False # No Odd number more than 1
f = factors5(x)
ps = powerset(f)
found = {y for y in {sum(i) for i in ps}
if 1 <= y < x}
return len(found) == x - 1
if __name__ == '__main__':
n = 333
p = [x for x in range(1, n + 1) if is_practical(x)]
print(f"There are {len(p)} Practical numbers from 1 to {n}:")
print(' ', str(p[:10])[1:-1], '...', str(p[-10:])[1:-1])
x = 666
print(f"\nSTRETCH GOAL: {x} is {'not ' if not is_practical(x) else }Practical.")</lang>
- Output:
There are 77 Practical numbers from 1 to 333: 1, 2, 4, 6, 8, 12, 16, 18, 20, 24 ... 288, 294, 300, 304, 306, 308, 312, 320, 324, 330 STRETCH GOAL: 666 is Practical.
Python: Alternate version
This version has an optimisation that might prove ``much`` faster in some cases but slower than the above in others.
A number with a large number of factors, f has 2**len(f) sets in its powerset. 672 for example has 23 factors and so 8_388_608 sets in its powerset.
Just taking the sets as they are generated and stopping when we first know that 672 is Practical needs just the first 28_826 or 0.34% of
the sets. 720, another Practical number needs just 0.01% of its half a billion sets to prove it is Practical.
Note however that if the number is ```not```
Composition of pure functions
Or as a composition of pure functions, including a more narrowly targeted test for subsets with a given sum:
<lang python>Practical numbers
from itertools import accumulate, chain, groupby, product from math import floor, sqrt from operator import mul from functools import reduce
- isPractical :: Int -> Bool
def isPractical(n):
True if n is a Practical number
(a member of OEIS A005153)
ds = properDivisors(n)
return all(map(
sumOfAnySubset(ds),
range(1, n)
))
- sumOfAnySubset :: [Int] -> Int -> Bool
def sumOfAnySubset(xs):
True if any subset of xs sums to n.
def go(n):
if n in xs:
return True
else:
total = sum(xs)
if n == total:
return True
elif n < total:
h, *t = reversed(xs)
d = n - h
return d in t or (
d > 0 and sumOfAnySubset(t)(d)
) or sumOfAnySubset(t)(n)
else:
return False
return go
- ------------------------- TEST -------------------------
- main :: IO ()
def main():
Practical numbers in the range [1..333],
and the OEIS A005153 membership of 666.
xs = [x for x in range(1, 334) if isPractical(x)]
print(
f'{len(xs)} OEIS A005153 numbers in [1..333]\n\n' + (
spacedTable(
chunksOf(10)([
str(x) for x in xs
])
)
)
)
print("\n")
for n in [666]:
print(
f'{n} is practical :: {isPractical(n)}'
)
- ----------------------- GENERIC ------------------------
- chunksOf :: Int -> [a] -> a
def chunksOf(n):
A series of lists of length n, subdividing the
contents of xs. Where the length of xs is not evenly
divible, the final list will be shorter than n.
def go(xs):
return [
xs[i:n + i] for i in range(0, len(xs), n)
] if 0 < n else None
return go
- primeFactors :: Int -> [Int]
def primeFactors(n):
A list of the prime factors of n.
def f(qr):
r = qr[1]
return step(r), 1 + r
def step(x):
return 1 + (x << 2) - ((x >> 1) << 1)
def go(x):
root = floor(sqrt(x))
def p(qr):
q = qr[0]
return root < q or 0 == (x % q)
q = until(p)(f)(
(2 if 0 == x % 2 else 3, 1)
)[0]
return [x] if q > root else [q] + go(x // q)
return go(n)
- properDivisors :: Int -> [Int]
def properDivisors(n):
The ordered divisors of n, excluding n itself.
def go(a, x):
return [a * b for a, b in product(
a,
accumulate(chain([1], x), mul)
)]
return sorted(
reduce(go, [
list(g) for _, g in groupby(primeFactors(n))
], [1])
)[:-1] if 1 < n else []
def listTranspose(xss):
Transposed matrix
def go(xss):
if xss:
h, *t = xss
return (
[[h[0]] + [xs[0] for xs in t if xs]] + (
go([h[1:]] + [xs[1:] for xs in t])
)
) if h and isinstance(h, list) else go(t)
else:
return []
return go(xss)
- until :: (a -> Bool) -> (a -> a) -> a -> a
def until(p):
The result of repeatedly applying f until p holds.
The initial seed value is x.
def go(f):
def g(x):
v = x
while not p(v):
v = f(v)
return v
return g
return go
- ---------------------- FORMATTING ----------------------
- spacedTable :: String -> String
def spacedTable(rows):
Tabulation with right-aligned cells
columnWidths = [
len(str(row[-1])) for row in listTranspose(rows)
]
def aligned(s, w):
return s.rjust(w, ' ')
return '\n'.join(
' '.join(
map(aligned, row, columnWidths)
) for row in rows
)
- MAIN ---
if __name__ == '__main__':
main()</lang>
- Output:
77 OEIS A005153 numbers in [1..333] 1 2 4 6 8 12 16 18 20 24 28 30 32 36 40 42 48 54 56 60 64 66 72 78 80 84 88 90 96 100 104 108 112 120 126 128 132 140 144 150 156 160 162 168 176 180 192 196 198 200 204 208 210 216 220 224 228 234 240 252 256 260 264 270 272 276 280 288 294 300 304 306 308 312 320 324 330 666 is practical :: True
Raku
<lang perl6>use Prime::Factor:ver<0.3.0+>;
sub is-practical ($n) {
my @proper-sums = $n.&proper-divisors.combinations».sum.unique.sort; +@proper-sums < $n-1 ?? False !! @proper-sums[^$n] eqv (^$n).list ?? True !! False
}
say "{+$_} matching numbers:\n{.batch(10)».fmt('%3d').join: "\n"}\n"
given [ (1..333).hyper(:32batch).grep: { is-practical($_) } ];
printf "%5s is practical? %s\n", $_, .&is-practical for 666, 6666, 66666;</lang>
- Output:
77 matching numbers: 1 2 4 6 8 12 16 18 20 24 28 30 32 36 40 42 48 54 56 60 64 66 72 78 80 84 88 90 96 100 104 108 112 120 126 128 132 140 144 150 156 160 162 168 176 180 192 196 198 200 204 208 210 216 220 224 228 234 240 252 256 260 264 270 272 276 280 288 294 300 304 306 308 312 320 324 330 666 is practical? True 6666 is practical? True 66666 is practical? False