%I #16 Jan 19 2024 04:54:25
%S 1,2,3,5,6,7,8,10,11,13,14,15,17,19,21,22,23,24,26,27,29,30,31,33,34,
%T 35,37,38,39,40,41,42,43,46,47,51,53,54,55,56,57,58,59,61,62,65,66,67,
%U 69,70,71,73,74,77,78,79,82,83,85,86,87,88,89,91,93,94,95
%N Numbers whose exponents in their prime factorization are either 1, 3, or both.
%C The asymptotic density of this sequence is zeta(6)/(zeta(2) * zeta(3)) * Product_{p prime} (1 + 2/p^3 - 1/p^4 + 1/p^5) = 0.68428692418686231814196872579121808347231273672316377728461822629005... (Cohen, 1962).
%C First differs from A036537 at n = 89. A036537(89) = 128 = 2^7 is not a term of this sequence.
%H Amiram Eldar, <a href="/A336591/b336591.txt">Table of n, a(n) for n = 1..10000</a>
%H Eckford Cohen, <a href="https://projecteuclid.org/euclid.pjm/1103036708">Arithmetical notes. III. Certain equally distributed sets of integers</a>, Pacific Journal of Mathematics, No. 12, Vol. 1 (1962), pp. 77-84.
%e 1 is a term since it has no exponents, and thus it has no exponent that is not 1 or 3.
%e 2 is a term since 2 = 2^1 has only the exponent 1 in its prime factorization.
%e 24 is a term since 24 = 2^3 * 3^1 has the exponents 1 and 3 in its prime factorization.
%t seqQ[n_] := AllTrue[FactorInteger[n][[;;,2]], MemberQ[{1, 3}, #] &]; Select[Range[100], seqQ]
%o (Python)
%o from itertools import count, islice
%o from sympy import factorint
%o def A336591_gen(startvalue=1): # generator of terms >= startvalue
%o return filter(lambda n:all(e==1 or e==3 for e in factorint(n).values()),count(max(startvalue,1)))
%o A336591_list = list(islice(A336591_gen(),20)) # _Chai Wah Wu_, Jun 22 2023
%Y Intersection of A046100 and A036537.
%Y Intersection of A046100 and A268335.
%Y A005117 and A062838 are subsequences.
%Y Cf. A068468.
%K nonn
%O 1,2
%A _Amiram Eldar_, Jul 26 2020