%I #33 Feb 12 2023 20:51:59
%S 1,8,9,18,20,27,32,36,42,44,45,49,50,54,63,68,72,78,80,81,84,90,92,99,
%T 105,108,110,114,116,117,125,126,128,135,144,153,156,162,164,168,169,
%U 170,171,176,180,186,188,189,195,198,200,207,210,212,216,222,225,228,230
%N Nonprime terms of A359059.
%C Subsequence of A359059 after we eliminate primes.
%e 8 is a term because 3|(4+2+12).
%e 9 is a term because 3|(6+3+12).
%e 18 is a term because 3|(6+8+36).
%e 20 is a term because 3|(8+10+36).
%e 27 is a term because 3|(18+3+36).
%p filter:= proc(n) local F,p,ph,r,ps;
%p F:= numtheory:-factorset(n);
%p if F = {n} then return false fi;
%p ph:= n * mul((p-1)/p, p = F);
%p r:= convert(F,`*`);
%p ps:= n * mul((p+1)/p, p = F);
%p (ph+r+ps) mod 3 = 0
%p end proc:
%p select(filter, [$1..1000]); # _Robert Israel_, Dec 20 2022
%t q[n_] := Module[{f = FactorInteger[n], p, e}, p = f[[;; , 1]]; e = f[[;; , 2]]; Divisible[Times @@ ((p - 1)*p^(e - 1)) + Times @@ p + Times @@ ((p + 1)*p^(e - 1)), 3]]; Select[Range[230], ! PrimeQ[#] && q[#] &] (* _Amiram Eldar_, Dec 20 2022 *)
%o (Python)
%o from sympy.ntheory.factor_ import totient
%o from sympy import isprime, primefactors, prod
%o def rad(n): return 1 if n < 2 else prod(primefactors(n))
%o def psi(n):
%o plist = primefactors(n)
%o return n*prod(p+1 for p in plist)//prod(plist)
%o # Output display terms.
%o for n in range(1,231):
%o if(False == isprime(n)):
%o if(0 == (totient(n) + rad(n) + psi(n)) % 3):
%o print(n, end = ", ")
%o (PARI) isok(m) = !isprime(m) && (((eulerphi(m) + factorback(factorint(m)[, 1]) + m*sumdiv(m, d, moebius(d)^2/d)) % 3) == 0); \\ _Michel Marcus_, Dec 27 2022
%Y Cf. A000010 (phi), A001615 (psi), A007947 (rad).
%Y Subsequence of A359059.
%K nonn
%O 1,2
%A _Torlach Rush_, Dec 15 2022