A363227 - OEIS

%I #60 Nov 12 2023 13:04:08

%S 1,2,3,4,6,8,9,10,12,14,15,16,18,20,21,24,27,28,30,32,36,40,42,44,45,

%T 48,50,52,54,56,60,63,64,66,70,72,75,78,80,81,84,88,90,96,98,99,100,

%U 102,104,105,108,110,112,114,117,120,126,128,130,132,135,136,138,140,144,147,150

%N Positive integers m such that every k with 1 <= k <= m is a linear combination of distinct divisors of m with coefficients +1 or -1.

%C All practical numbers are terms of this sequence.

%C Referred to as "semi-practical" numbers in the comments of the linked Numberphile video.

%H Alois P. Heinz, <a href="/A363227/b363227.txt">Table of n, a(n) for n = 1..10000</a>

%H James Grime, <a href="https://www.youtube.com/watch?v=IlZOLwf87gM">Practical Numbers</a>, Numberphile video (2023).

%p q:= proc(n) uses numtheory, ListTools; local b, l, s;

%p       l:= [divisors(n)[]]; s:= PartialSums(l);

%p       b:= proc(m, i) option remember; m=0 or i>0 and m<=s[i] and

%p            (b(m, i-1) or b(abs(m-l[i]), i-1) or b(m+l[i], i-1))

%p           end;

%p       andmap(x-> b(x, nops(l)), [$1..n])

%p     end:

%p select(q, [$1..256])[];  # _Alois P. Heinz_, May 23 2023

%t q[n_] := Module[{b, l, s}, l = Divisors[n]; s = Accumulate[l]; b[m_, i_] := b[m, i] = m == 0 || i > 0 && m <= s[[i]] && (b[m, i-1] || b[Abs[m - l[[i]]], i-1] || b[m+l[[i]], i-1]); AllTrue[Range[n], b[#, Length[l]]&]];

%t Select[Range[256], q] (* _Jean-François Alcover_, Nov 12 2023, after _Alois P. Heinz_ *)

%o (Python)

%o import math

%o def divisorGenerator(n):

%o   large_divisors = []

%o   for i in range(1, int(math.sqrt(n) + 1)):

%o     if n % i == 0:

%o       yield i

%o       if i*i != n:

%o         large_divisors.append(int(n / i))

%o     for divisor in reversed(large_divisors):

%o       yield divisor

%o from itertools import chain, combinations

%o def all_combinations(iterable,n):

%o   s = list(iterable)

%o   for sumset in chain.from_iterable(combinations(s, r) for r in range(len(s)+1)):

%o     remaining = list(set(s).symmetric_difference(set(list(sumset))))

%o     for subtractset in chain.from_iterable(combinations(remaining, r) for r in range(len(remaining)+1)):

%o       value = sum(list(sumset))-sum(list(subtractset))

%o       if value>0 and value<=n:

%o         yield value

%o def is_A363227(n):

%o   return len(set(all_combinations(divisorGenerator(n),n)))==n

%o max_n = 250

%o print([x for x in range(max_n+1) if is_A363227(x)])

%o (Python)

%o from itertools import count, islice

%o from sympy import divisors

%o def A363227_gen(startvalue=1): # generator of terms >= startvalue

%o     for m in count(max(startvalue,1)):

%o         c = {0}

%o         for d in divisors(m,generator=True):

%o             c |= {a+d for a in c}|{a-d for a in c}

%o         if all(k in c for k in range(1,m+1)):

%o             yield m

%o A363227_list = list(islice(A363227_gen(),20)) # _Chai Wah Wu_, Jul 04 2023

%Y Cf. A005153 (practical numbers).

%K nonn

%O 1,2

%A _Guy Ziv_, May 21 2023