login
Nonprime terms of A359059.
0

%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