A339509 - OEIS

%I #21 Jul 09 2021 10:21:15

%S 1,2,3,4,6,7,9,10,14,18,20,21,29,30,32,36,44,45,61,62,70,74,76,77,109,

%T 125,127,191,199,200,232,233,249,253,255,271,399,400,402,406,470,471,

%U 503,504,512,640,642,643,899,963,1219,1223,1231,1232,1744,1760,1888

%N Number of subsets of {1..n} whose elements have the same greatest prime factor.

%H Alois P. Heinz, <a href="/A339509/b339509.txt">Table of n, a(n) for n = 0..10000</a>

%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/GreatestPrimeFactor.html">Greatest Prime Factor</a>

%e a(8) = 14 subsets: {}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {2, 4}, {2, 8}, {3, 6}, {4, 8} and {2, 4, 8}.

%p b:= proc(n) option remember; `if`(n<2, 0,

%p       b(n-1)+x^max(numtheory[factorset](n)))

%p     end:

%p a:= n-> `if`(n<2, n+1, (p-> 2+add(2^

%p     coeff(p, x, i)-1, i=2..degree(p)))(b(n))):

%p seq(a(n), n=0..70);  # _Alois P. Heinz_, Dec 07 2020

%t b[n_] := b[n] = If[n < 2, 0,

%t      b[n - 1] + x^Max[FactorInteger[n][[All, 1]]]];

%t a[n_] := If[n < 2, n + 1, Function[p, 2 + Sum[2^

%t      Coefficient[p, x, i] - 1, {i, 2, Exponent[p, x]}]][b[n]]];

%t Table[a[n], {n, 0, 70}] (* _Jean-François Alcover_, Jul 09 2021, after _Alois P. Heinz_ *)

%o (Python)

%o from sympy import primefactors

%o def test(n):

%o     if n<2: return n

%o     return max(primefactors(n))

%o def a(n):

%o     tests = [test(i) for i in range(n+1)]

%o     return sum(2**tests.count(v)-1 for v in set(tests))

%o print([a(n) for n in range(57)]) # _Michael S. Branicky_, Dec 11 2020

%Y Cf. A006530, A339510.

%K nonn

%O 0,2

%A _Ilya Gutkovskiy_, Dec 07 2020

%E a(24)-a(56) from _Alois P. Heinz_, Dec 07 2020