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
APL
<lang APL>pract ← ∧/⍳∊(+⌿⊢×[1](2/⍨≢)⊤(⍳2*≢))∘(⍸0=⍳|⊢)</lang>
- Output:
<lang APL> ⍸pract¨⍳333 ⍝ Which numbers from 1 to 333 are practical? 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
pract 666 ⍝ Is 666 practical?
1</lang>
C#
<lang csharp>using System.Collections.Generic; using System.Linq; using static System.Console;
class Program {
static bool soas(int n, IEnumerable<int> f) {
if (n <= 0) return false; if (f.Contains(n)) return true;
switch(n.CompareTo(f.Sum())) { case 1: return false; case 0: return true;
case -1: var rf = f.Reverse().ToList(); var d = n - rf[0]; rf.RemoveAt(0);
return soas(d, rf) || soas(n, rf); } return true; }
static bool ip(int n) { var f = Enumerable.Range(1, n >> 1).Where(d => n % d == 0).ToList();
return Enumerable.Range(1, n - 1).ToList().TrueForAll(i => soas(i, f)); }
static void Main() {
int c = 0, m = 333; for (int i = 1; i <= m; i += i == 1 ? 1 : 2)
if (ip(i) || i == 1) Write("{0,3} {1}", i, ++c % 10 == 0 ? "\n" : "");
Write("\nFound {0} practical numbers between 1 and {1} inclusive.\n", c, m);
do Write("\n{0,5} is a{1}practical number.",
m = m < 500 ? m << 1 : m * 10 + 6, ip(m) ? " " : "n im"); while (m < 1e4); } }</lang>
- Output:
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 Found 77 practical numbers between 1 and 333 inclusive. 666 is a practical number. 6666 is a practical number. 66666 is an impractical number.
C++
<lang cpp>#include <algorithm>
- include <iostream>
- include <numeric>
- include <sstream>
- include <vector>
// Returns true if any subset of f sums to n. bool sum_of_any_subset(int n, const std::vector<int>& f) {
if (f.empty())
return false;
if (std::find(f.begin(), f.end(), n) != f.end())
return true;
int total = std::accumulate(f.begin(), f.end(), 0);
if (n == total)
return true;
if (n > total)
return false;
int d = n - f.back();
std::vector<int> g(f.begin(), f.end() - 1);
return std::find(g.begin(), g.end(), d) != g.end() ||
(d > 0 && sum_of_any_subset(d, g)) || sum_of_any_subset(n, g);
}
// Returns the proper divisors of n. std::vector<int> factors(int n) {
std::vector<int> f{1};
for (int i = 2; i * i <= n; ++i) {
if (n % i == 0) {
f.push_back(i);
if (i * i != n)
f.push_back(n / i);
}
}
std::sort(f.begin(), f.end());
return f;
}
bool is_practical(int n) {
std::vector<int> f = factors(n);
for (int i = 1; i < n; ++i) {
if (!sum_of_any_subset(i, f))
return false;
}
return true;
}
std::string shorten(const std::vector<int>& v, size_t n) {
std::ostringstream out;
size_t size = v.size(), i = 0;
if (n > 0 && size > 0)
out << v[i++];
for (; i < n && i < size; ++i)
out << ", " << v[i];
if (size > i + n) {
out << ", ...";
i = size - n;
}
for (; i < size; ++i)
out << ", " << v[i];
return out.str();
}
int main() {
std::vector<int> practical;
for (int n = 1; n <= 333; ++n) {
if (is_practical(n))
practical.push_back(n);
}
std::cout << "Found " << practical.size() << " practical numbers:\n"
<< shorten(practical, 10) << '\n';
std::cout << std::boolalpha;
for (int n : {666, 6666, 66666, 672, 720})
std::cout << "is_practical(" << n << "): " << is_practical(n) << '\n';
return 0;
}</lang>
- 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
Delphi
See Pascal.
J
Borrowed from the Proper divisors#J page:
<lang J>factors=: [: /:~@, */&>@{@((^ i.@>:)&.>/)@q:~&__ properDivisors=: factors -. ]</lang>
Borrowed from the Power set#J page:
<lang J>ps=: #~ 2 #:@i.@^ #</lang>
Implementation:
<lang J>isPrac=: ( -:&# i. -. 0,+/"1@(ps ::empty)@properDivisors)"0</lang>
Task examples:
<lang J> +/ isPrac 1+i.333 NB. count practical numbers 77
(#~ isPrac) 1+i.333 NB. list them
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 30...
isPrac 666 NB. test
1</lang>
Julia
<lang julia>using Primes
""" proper divisors of n """ function proper_divisors(n)
f = [one(n)]
for (p,e) in factor(n)
f = reduce(vcat, [f*p^j for j in 1:e], init=f)
end
pop!(f)
return f
end
""" return true if any subset of f sums to n. """ function sumofanysubset(n, f)
n in f && return true total = sum(f) n == total && return true n > total && return false rf = reverse(f) d = n - popfirst!(rf) return (d in rf) || (d > 0 && sumofanysubset(d, rf)) || sumofanysubset(n, rf)
end
function ispractical(n)
n == 1 && return true isodd(n) && return false f = proper_divisors(n) return all(i -> sumofanysubset(i, f), 1:n-1)
end
const prac333 = filter(ispractical, 1:333) println("There are ", length(prac333), " practical numbers between 1 to 333 inclusive.") println("Start and end: ", filter(x -> x < 25, prac333), " ... ", filter(x -> x > 287, prac333), "\n") for n in [666, 6666, 66666, 222222]
println("$n is ", ispractical(n) ? "" : "not ", "a practical number.")
end
</lang>
- Output:
There are 77 practical numbers between 1 to 333 inclusive. Start and end: [1, 2, 4, 6, 8, 12, 16, 18, 20, 24] ... [288, 294, 300, 304, 306, 308, 312, 320, 324, 330] 666 is a practical number. 6666 is a practical number. 66666 is not a practical number. 222222 is a practical number.
Pascal
simple brute force.Marking sum of divs by shifting the former sum by the the next divisor.
SumAllSetsForPractical tries to break as soon as possible.Should try to check versus https://en.wikipedia.org/wiki/Practical_number#Characterization_of_practical_numbers
...σ denotes the sum of the divisors of x. For example, 2 × 3^2 × 29 × 823 = 429606 is practical, because the inequality above holds for each of its prime factors: 3 ≤ σ(2) + 1 = 4, 29 ≤ σ(2 × 3^2) + 1 = 40, and 823 ≤ σ(2 × 3^2 × 29) + 1 = 1171.
<lang pascal>program practicalnumbers; {$IFDEF FPC}
{$MODE DELPHI}{$OPTIMIZATION ON,ALL}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils
{$IFNDEF FPC}
,Windows
{$ENDIF}
;
const
LOW_DIVS = 0; HIGH_DIVS = 2048 - 1;
type
tdivs = record DivsVal: array[LOW_DIVS..HIGH_DIVS] of Uint32; DivsMaxIdx, DivsNum, DivsSumProp: NativeUInt; end;
var
Divs: tDivs; HasSum: array of byte;
procedure GetDivisors(var Divs: tdivs; n: Uint32); //calc all divisors,keep sorted var
i, quot, ug, og: UInt32; sum: UInt64;
begin
with Divs do begin DivsNum := n; sum := 0; ug := 0; og := HIGH_DIVS; i := 1;
while i * i < n do
begin
quot := n div i;
if n - quot * i = 0 then
begin
DivsVal[og] := quot;
Divs.DivsVal[ug] := i;
inc(sum, quot + i);
dec(og);
inc(ug);
end;
inc(i);
end;
if i * i = n then
begin
DivsVal[og] := i;
inc(sum, i);
dec(og);
end;
//move higher divisors down
while og < high_DIVS do
begin
inc(og);
DivsVal[ug] := DivsVal[og];
inc(ug);
end;
DivsMaxIdx := ug - 2;
DivsSumProp := sum - n;
end; //with
end;
function SumAllSetsForPractical(Limit: Uint32): boolean; //mark sum and than shift by next divisor == add //for practical numbers every sum must be marked var
hs0, hs1: pByte; idx, j, loLimit, maxlimit, delta: NativeUint;
begin
Limit := trunc(Limit * (Limit / Divs.DivsSumProp));
loLimit := 0;
maxlimit := 0;
hs0 := @HasSum[0];
hs0[0] := 1; //empty set
for idx := 0 to Divs.DivsMaxIdx do
begin
delta := Divs.DivsVal[idx];
hs1 := @hs0[delta];
for j := maxlimit downto 0 do
hs1[j] := hs1[j] or hs0[j];
maxlimit := maxlimit + delta;
while hs0[loLimit] <> 0 do
inc(loLimit);
//IF there is a 0 < delta, it will never be set
//IF there are more than the Limit set,
//it will be copied by the following Delta's
if (loLimit < delta) or (loLimit > Limit) then
Break;
end;
result := (loLimit > Limit);
end;
function isPractical(n: Uint32): boolean; var
i: NativeInt; sum: NativeUInt;
begin
if n = 1 then EXIT(True); if ODD(n) then EXIT(false); if (n > 2) and not ((n mod 4 = 0) or (n mod 6 = 0)) then EXIT(false);
GetDivisors(Divs, n);
i := n - 1;
sum := Divs.DivsSumProp;
if sum < i then
result := false
else
begin
if length(HasSum) > sum + 1 + 1 then
FillChar(HasSum[0], sum + 1, #0)
else
begin
setlength(HasSum, 0);
setlength(HasSum, sum + 8 + 1);
end;
result := SumAllSetsForPractical(i);
end;
end;
procedure OutIsPractical(n: nativeInt); begin
if isPractical(n) then writeln(n, ' is practical') else writeln(n, ' is not practical');
end;
const
ColCnt = 10; MAX = 333;
var
T0: Int64; n, col, count: NativeInt;
begin
col := ColCnt;
count := 0;
for n := 1 to MAX do
if isPractical(n) then
begin
write(n: 5);
inc(count);
dec(col);
if col = 0 then
begin
writeln;
col := ColCnt;
end;
end;
writeln;
writeln('There are ', count, ' pratical numbers from 1 to ', MAX);
writeln;
T0 := GetTickCount64;
OutIsPractical(666);
OutIsPractical(6666);
OutIsPractical(66666);
OutIsPractical(954432);
OutIsPractical(720);
OutIsPractical(5384);
OutIsPractical(1441440);
writeln(Divs.DivsNum, ' has ', Divs.DivsMaxIdx + 1, ' proper divisors');
writeln((GetTickCount64 - T0) / 1000: 10: 3, ' s');
T0 := GetTickCount64;
OutIsPractical(99998640);
writeln(Divs.DivsNum, ' has ', Divs.DivsMaxIdx + 1, ' proper divisors ');
writeln((GetTickCount64 - T0) / 1000: 10: 3, ' s');
T0 := GetTickCount64;
OutIsPractical(99998640);
writeln(Divs.DivsNum, ' has ', Divs.DivsMaxIdx + 1, ' proper divisors ');
writeln((GetTickCount64 - T0) / 1000: 10: 3, ' s');
setlength(HasSum, 0);
{$IFNDEF UNIX} readln; {$ENDIF}
end. </lang>
- Output:
TIO.RUN.
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
There are 77 pratical numbers from 1 to 333
666 is practical
6666 is practical
66666 is not practical
954432 is not practical
720 is practical
5384 is not practical
1441440 is practical
1441440 has 287 proper divisors
0.017 s
99998640 is not practical
99998640 has 119 proper divisors
0.200 s // with reserving memory
99998640 is not practical
99998640 has 119 proper divisors
0.081 s // already reserved memory
Real time: 0.506 s CPU share: 87.94 %
alternative
Now without generating sum of allset. <lang pascal>program practicalnumbers2;
{$IFDEF FPC}
{$MODE DELPHI}{$OPTIMIZATION ON,ALL}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} uses
SysUtils;
type
tdivs = record DivsVal: array[0..1024 - 1] of Uint32; end;
var
Divs: tDivs;
function CheckIsPractical(var Divs: tdivs; n: Uint32): boolean; //calc all divisors,calc sum of divisors var sum: UInt64; i :NativeInt; quot,ug,sq,de: UInt32;
begin
with Divs do
begin
sum := 1;
ug := Low(tdivs.DivsVal);
i := 2;
sq := 4;
de := 5;
while sq < n do
begin
quot := n div i;
if n - quot * i = 0 then
begin
if sum + 1 < i then
EXIT(false);
Inc(sum, i);
DivsVal[ug] := quot;
Inc(ug);
end;
Inc(i);
sq += de;
de := de+2;
end;
if sq = n then
begin
if sum + 1 < i then
EXIT(false);
DivsVal[ug] := i;
Inc(sum, i);
Inc(ug);
end;
//check higher
while ug > 0 do
begin
Dec(ug);
i := DivsVal[ug];
if sum + 1 < i then
EXIT(false);
Inc(sum, i);
if sum >= n - 1 then
break;
end;
end;//with
result := true;
end;
function isPractical(n: Uint32): boolean;
begin
if n in [1,2] then
EXIT(True);
if ODD(n) then
EXIT(False);
Result := CheckIsPractical(Divs, n);
end;
procedure OutIsPractical(n: nativeInt);
begin
if isPractical(n) then
writeln(n, ' is practical')
else
writeln(n, ' is not practical');
end;
const
ColCnt = 10; MAX = 333;
var
T0 : int64; n, col, Count: NativeInt;
begin
col := ColCnt;
Count := 0;
for n := 1 to MAX do
if isPractical(n) then
begin
Write(n: 5);
Inc(Count);
Dec(col);
if col = 0 then
begin
writeln;
col := ColCnt;
end;
end;
writeln;
writeln('There are ', Count, ' pratical numbers from 1 to ', MAX);
writeln;
OutIsPractical(666); OutIsPractical(6666); OutIsPractical(66666); OutIsPractical(954432); OutIsPractical(720); OutIsPractical(5184); OutIsPractical(1441440); OutIsPractical(99998640);
T0 := GetTickCOunt64;
count := 0;
For n := 1 to 1000*1000 do
inc(count,Ord(isPractical(n)));
writeln('Count of practical numbers til 1,000,000 ',count,(GetTickCount64-t0)/1000:8:4,' s');
{$IFDEF WINDOWS}
readln;
{$ENDIF}
end. </lang>
- Output:
TIO.RUN
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
There are 77 pratical numbers from 1 to 333
666 is practical
6666 is practical
66666 is not practical
954432 is not practical
720 is practical
5184 is practical
1441440 is practical
99998640 is not practical
Count of practical numbers til 1,000,000 97385 2.1380 s
Real time: 2.277 s CPU share: 99.55 %
Perl
<lang perl>use strict; use warnings; use feature 'say'; use enum <False True>; use ntheory <divisors vecextract>; use List::AllUtils <sum uniq firstidx>;
sub proper_divisors {
return 1 if 0 == (my $n = shift); my @d = divisors($n); pop @d; @d
}
sub powerset_sums { uniq map { sum vecextract(\@_,$_) } 1..2**@_-1 }
sub is_practical {
my($n) = @_;
return True if $n == 1;
return False if 0 != $n % 2;
($n-2 == firstidx { $_ == $n-1 } powerset_sums(proper_divisors($n)) ) ? True : False;
}
my @pn; is_practical($_) and push @pn, $_ for 1..333; say @pn . " matching numbers:\n" . (sprintf "@{['%4d' x @pn]}", @pn) =~ s/(.{40})/$1\n/gr; say ; printf "%6d is practical? %s\n", $_, is_practical($_) ? 'True' : 'False' 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
Phix
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
- %% 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: Faster version
This version has an optimisation that proves much faster when testing a range of numbers for Practicality.
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.
The inner loop is sensitive to the order of factors passed to the powerset generator and experimentation shows that reverse sorting the factors saves the most computation.
An extra check on the sum of all factors has a minor positive effect too.
<lang python>def is_practical5(x: int) -> bool:
"""Practical number test with factor reverse sort and short-circuiting."""
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 = sorted(factors5(x), reverse=True)
if sum(f) < x - 1:
return False # Never get x-1
ps = powerset(f)
found = set()
for nps in ps:
if len(found) < x - 1:
y = sum(nps)
if 1 <= y < x:
found.add(y)
else:
break # Short-circuiting the loop.
return len(found) == x - 1
if __name__ == '__main__':
n = 333
print("\n" + is_practical5.__doc__)
p = [x for x in range(1, n + 1) if is_practical5(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.")
x = 5184
print(f"\nEXTRA GOAL: {x} is {'not ' if not is_practical(x) else }Practical.")</lang>
- Output:
Using the definition of factors5 from the simple case above then the following results are obtained.
Practical number test with factor reverse sort and short-circuiting. 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. EXTRA GOAL: 5184 is Practical.
5184, which is practical, has 34 factors!
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.
Around 1/8'th of the integers up to the 10_000'th Practical number have more than 22 factors but are not practical numbers themselves. (Some of these have 31 factors whilst being divisible by four or six so would need the long loop to complete)!
Composition of pure functions
<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 Callable, List
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[[int], bool]:
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 -------------------------
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)}'
)
- ----------------------- GENERIC ------------------------
def chunksOf(n: int) -> Callable[[List[str]], List[List[str]]]:
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[str]]) -> List[List[str]]:
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: Callable[[int], bool]) -> Callable[[int], bool]:
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 ----------------------
def spacedTable(rows: List[List[str]]) -> 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
)
- 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) {
return True if $n == 1;
return False if $n % 2;
my @proper = $n.&proper-divisors :sort;
return True if all( @proper.rotor(2 => -1).map: { .[0] / .[1] >= .5 } );
my @proper-sums = @proper.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, 672, 720;</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 672 is practical? True 720 is practical? True
Visual Basic .NET
<lang vbnet>Imports System.Collections.Generic, System.Linq, System.Console
Module Module1
Function soas(ByVal n As Integer, ByVal f As IEnumerable(Of Integer)) As Boolean
If n <= 0 Then Return False Else If f.Contains(n) Then Return True
Select Case n.CompareTo(f.Sum())
Case 1 : Return False : Case 0 : Return True
Case -1 : Dim rf As List(Of Integer) = f.Reverse().ToList() : Dim D as Integer = n - rf(0)
rf.RemoveAt(0) : Return soas(d, rf) OrElse soas(n, rf)
End Select : Return true
End Function
Function ip(ByVal n As Integer) As Boolean
Dim f As IEnumerable(Of Integer) = Enumerable.Range(1, n >> 1).Where(Function(d) n Mod d = 0).ToList()
Return Enumerable.Range(1, n - 1).ToList().TrueForAll(Function(i) soas(i, f))
End Function
Sub Main()
Dim c As Integer = 0, m As Integer = 333, i As Integer = 1 : While i <= m
If ip(i) OrElse i = 1 Then c += 1 : Write("{0,3} {1}", i, If(c Mod 10 = 0, vbLf, ""))
i += If(i = 1, 1, 2) : End While
Write(vbLf & "Found {0} practical numbers between 1 and {1} inclusive." & vbLf, c, m)
Do : m = If(m < 500, m << 1, m * 10 + 6)
Write(vbLf & "{0,5} is a{1}practical number.", m, If(ip(m), " ", "n im")) : Loop While m < 1e4
End Sub
End Module</lang>
- Output:
Same as C#
Wren
<lang ecmascript>import "/math" for Int, Nums
var powerset // recursive powerset = Fn.new { |set|
if (set.count == 0) return [[]]
var head = set[0]
var tail = set[1..-1]
return powerset.call(tail) + powerset.call(tail).map { |s| [head] + s }
}
var isPractical = Fn.new { |n|
if (n == 1) return true
var divs = Int.properDivisors(n)
var subsets = powerset.call(divs)
var found = List.filled(n, false)
var count = 0
for (subset in subsets) {
var sum = Nums.sum(subset)
if (sum > 0 && sum < n && !found[sum]) {
found[sum] = true
count = count + 1
if (count == n - 1) return true
}
}
return false
}
System.print("In the range 1..333, there are:") var count = 0 var practical = [] for (i in 1..333) {
if (isPractical.call(i)) {
count = count + 1
practical.add(i)
}
} System.print(" %(count) practical numbers") System.print(" The first ten are %(practical[0..9])") System.print(" The final ten are %(practical[-10..-1])") System.print("\n666 is practical: %(isPractical.call(666))")</lang>
- Output:
In the range 1..333, there are: 77 practical numbers The first ten are [1, 2, 4, 6, 8, 12, 16, 18, 20, 24] The final ten are [288, 294, 300, 304, 306, 308, 312, 320, 324, 330] 666 is practical: true