%I #15 Sep 15 2024 02:42:45
%S 1,64,243,441,1764,9800,15552,28224,41616,60516,82369,88200,189728,
%T 226576,329476,336200,648675,741321,968256,1317904,1428025,1707552,
%U 1943236,2039184,2056356,2381400,2446227,2798929,2965284,2986568,4372281,5189400,5271616,6508832
%N Powerful numbers (A001694) whose sum of powerful divisors (including 1) is also powerful.
%C Numbers k such that A112526(k) = A112526(A183097(k)) = 1.
%H Amiram Eldar, <a href="/A349109/b349109.txt">Table of n, a(n) for n = 1..12154</a> (terms below 10^19)
%e 64 = 2^6 is a term since it is powerful and the sum of its powerful divisors, A183097(64) = 1 + 4 + 8 + 16 + 32 + 64 = 125 = 5^3 is also powerful.
%t powQ[n_] := n == 1 || AllTrue[FactorInteger[n][[;;,2]], # > 1 &]; f[p_, e_] := (p^(e + 1) - 1)/(p - 1) - p; s[1] = 1; s[n_] := Times @@ f @@@ FactorInteger[n]; q[n_] := powQ[n] && powQ[s[n]]; Select[Range[7*10^6], q]
%o (PARI) isok(n) = ispowerful(n) && ispowerful(sumdiv(n, d, d*ispowerful(d))); \\ _Michel Marcus_, Nov 08 2021
%o (PARI) is(k) = {my(f = factor(k)); ispowerful(f) && ispowerful(prod(i = 1, #f~, (f[i,1]^(f[i,2]+1) - 1)/(f[i,1] - 1) - f[i,1]));} \\ _Amiram Eldar_, Sep 14 2024
%Y Cf. A001694, A112526, A180090, A183097, A337044, A337045, A349110.
%K nonn,changed
%O 1,2
%A _Amiram Eldar_, Nov 08 2021