%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