Practical numbers

Practical numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

Translation of: Python – (the composition of functions version)

function sum_of_any_subset(integer n, sequence f)
    -- return true if any subset of f sums to n.
    if find(n,f) then return true end if
    integer total = sum(f)
    if n=total then return true
    elsif n>total then return false end if
    integer d = n-f[$]
    f = f[1..$-1]
    return find(d,f)
        or (d>0 and sum_of_any_subset(d, f))
        or sum_of_any_subset(n, f)
end function

function is_practical(integer n)
    sequence f = factors(n,-1)
    for i=1 to n-1 do
        if not sum_of_any_subset(i,f) then return false end if
    end for
    return true
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\n",{n,is_practical(n)})
end procedure
papply({666,6666,66666,672,720},stretch)

Output:

Found 77 practical numbers:
  1,   2,   4,   6,   8,  12,  16,  18,  20,  24, ..., 288, 294, 300, 304, 306, 308, 312, 320, 324, 330

is_practical(666):true
is_practical(6666):true
is_practical(66666):false
is_practical(672):true
is_practical(720):true

Python

Python: Straight forward implementation

<lang python>from itertools import chain, cycle, accumulate, combinations from typing import List, Tuple

1. %% 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]

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))

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 ultimately Practical, with a large number of factors, then the loop to find this is slower.

<lang python>def is_practical2(x: int, __switch=23) -> 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.

   Can short-circuit summations for x with large number of factors >= __switch
   """
   if x == 1:
       return True
   if x % 2:
       return False  # No Odd number more than 1
   mult_4_or_6 = (x % 4 == 0) or (x % 6 == 0)
   if x > 2 and not mult_4_or_6:
       return False  # If > 2 then must be a divisor of 4 or 6
   f = factors5(x)
   ps = powerset(f)    # = 2**len(f) items

   if len(f) < __switch:
       found = {y for y in {sum(i) for i in ps}
                if 1 <= y < x}
   else:
       found = set()
       while len(found) < x - 1:
           y = sum(next(ps))
           if 1 <= y < x:
               found.add(y)

   return len(found) == x - 1</lang>

Output:

If the above function replaces that in the simple case above then the same results are produced.

672, which is practical and has 23 factors computes 200 times faster.

A little further investigation shows that you need to get to 3850, for the first example of a number with 23 or more factors that is not Practical.

Composition of pure functions

A composition of pure functions, including a more narrowly targeted test for subsets with a given sum in the ranges covered here.

Note that no type hints are added for the compiler in this example. The additional work of adding such hints can pay off with larger projects, particularly those involving several collaborators, but the cost/benefit (work for the coder, run-time startup cost) is known to be rather less clear at a small scale.

Clarity about the return type semantics, is however, very useful when reasoning about pure functions, and my personal approach is to add light informal comments about the type, in a Hindley Milner idiom, which lends itself well to brief and clean notes on the type of curried functions, which are more easily composable, especially with higher order functions, and which I generally prefer to use.

I've been asked (always by the same person :-) why I don't find the idiom of Python compiler type-hints a good match for my semantic type comments, and the answer is essentially that the compiler type hints are not a clear or helpful notation for this purpose – not just because they generally involve more typing and visual noise, but also, and in particular, because with curried functions the compiler hint notation becomes swamped by use of the cognitively redundant `Callable` keyword, which degrades clarity, and imposes burden, for the human reader.

<lang python>Practical numbers

from itertools import accumulate, chain, groupby, product from math import floor, sqrt from operator import mul from functools import reduce

1. 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)
   ))

1. 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

1. ------------------------- TEST -------------------------
2. 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)}'
       )

1. ----------------------- GENERIC ------------------------

1. 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

1. 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)

1. 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 []

1. listTranspose :: a -> a

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)

1. 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

1. ---------------------- FORMATTING ----------------------

1. 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
   )

1. 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

Semantic type comments replaced by compiler type hints

Note the the prominence of the Callable keyword, and the loss of specificity for the reader, in the uses of the Any type symbol.

Semantic clarity for the human reader could be improved by restoring the informal – but cleaner, and more informative – type comments which this version deletes from the original version above. <lang python>Practical numbers

from itertools import accumulate, chain, groupby, product from math import floor, sqrt from operator import mul from functools import reduce from typing import Any, List, Callable

def isPractical(n: int) -> bool:

   True if n is a Practical number
      (a member of OEIS A005153)
   
   ds = properDivisors(n)
   return all(map(
       sumOfAnySubset(ds),
       range(1, n)
   ))

def sumOfAnySubset(xs: List[int]) -> Callable[[Any], Any]:

   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

1. ------------------------- TEST -------------------------

def main() -> None:

   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)}'
       )

1. ----------------------- GENERIC ------------------------

def chunksOf(n: int) -> Callable[[Any], Any]:

   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

def primeFactors(n: int) -> List[int]:

   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)

def properDivisors(n: int) -> List[int]:

   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: List[List[int]]) -> List[List[int]]:

   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)

def until(p:bool) -> Any:

   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

1. ---------------------- FORMATTING ----------------------

def spacedTable(rows: List[List[int]]) -> str:

   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
   )

1. 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