%I #25 Feb 25 2021 08:31:32
%S 4,8,12,16,20,24,25,27,28,32,36,40,44,48,49,50,52,54,56,60,64,68,72,
%T 75,76,80,81,84,88,92,96,98,100,104,108,112,116,120,121,124,125,128,
%U 132,135,136,140,144,147,148,150,152,156,160,162,164,168,169,172,175,176
%N Positive numbers divisible by 3^3 or by the square of some other prime.
%C Or numbers n>=4 having a divisor k^2>=4 such that n and n/k^2 equal modulo 3.
%C All positive multiples of 4 are in the sequence.
%C Or numbers n such that there is a smaller positive number j == n (mod 3) such that sqrt(j*n) is integer. The smallest such j corresponds to the greatest k; or, the same, j = 3*A007913(n/3), if n is divisible by 3 and otherwise j=A007913(n).
%C Or complement to the sequence: S, 3*S and 9*S, where S denotes the sequence of the squarefree numbers not divisible by 3.
%H Amiram Eldar, <a href="/A274141/b274141.txt">Table of n, a(n) for n = 1..10000</a>
%F Let A(x) be the number of a(n)<=x. Then A(x)~(1 - 6.5/Pi^2)*x = 0.34141230...*x as x goes to infinity.
%t Select[Range[200], (e = IntegerExponent[#, 3]) > 2 || ! SquareFreeQ[#/3^e] &] (* _Amiram Eldar_, Feb 25 2021 *)
%o (PARI) isok(n) = (((v=valuation(n, 3)) >= 3) || (((m = n/3^v) > 1) && (vecmax((factor(m))[,2]) >=2))); \\ _Michel Marcus_, Jun 12 2016
%Y Cf. A046790, A046791, A007913.
%K nonn
%O 1,1
%A _Vladimir Shevelev_, Jun 11 2016