A359062 Nonprime terms of A359059.

1, 8, 9, 18, 20, 27, 32, 36, 42, 44, 45, 49, 50, 54, 63, 68, 72, 78, 80, 81, 84, 90, 92, 99, 105, 108, 110, 114, 116, 117, 125, 126, 128, 135, 144, 153, 156, 162, 164, 168, 169, 170, 171, 176, 180, 186, 188, 189, 195, 198, 200, 207, 210, 212, 216, 222, 225, 228, 230

Offset: 1

Author

Torlach Rush, Dec 15 2022

Keywords

• nonn

Comments

Subsequence of A359059 after we eliminate primes.

Examples

			8 is a term because 3|(4+2+12).
9 is a term because 3|(6+3+12).
18 is a term because 3|(6+8+36).
20 is a term because 3|(8+10+36).
27 is a term because 3|(18+3+36).
		

Crossrefs

Cf. A000010 (phi), A001615 (psi), A007947 (rad).

Subsequence of A359059.

Programs

Maple

filter:= proc(n) local F,p,ph,r,ps;
    F:= numtheory:-factorset(n);
      if F = {n} then return false fi;
      ph:= n * mul((p-1)/p, p = F);
    r:= convert(F,`*`);
    ps:= n * mul((p+1)/p, p = F);
    (ph+r+ps) mod 3 = 0
end proc:
select(filter, [$1..1000]); # Robert Israel, Dec 20 2022

Mathematica

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 *)

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

Python

from sympy.ntheory.factor_ import totient
from sympy import isprime, primefactors, prod
def rad(n): return 1 if n < 2 else prod(primefactors(n))
def psi(n):
    plist = primefactors(n)
    return n*prod(p+1 for p in plist)//prod(plist)
# Output display terms.
for n in range(1,231):
    if(False == isprime(n)):
        if(0 == (totient(n) + rad(n) + psi(n)) % 3):
            print(n, end = ", ")