cp's OEIS Frontend

A331701 Prime powers (A025475) that can be represented as a sum of two prime powers.

%I A331701 #13 Feb 03 2020 16:45:41
%S A331701 8,9,16,25,32,64,81,125,128,256,512,1024,2048,4096,5041,8192,16384,
%T A331701 32768,65536,131072,262144,524288,1048576,2097152,4194304,8388608,
%U A331701 16777216,33554432,67108864,134217728,268435456,536870912,1073741824,2147483648,4294967296,8589934592
%N A331701 Prime powers (A025475) that can be represented as a sum of two prime powers.
%C A331701 A000079 is a subsequence, starting from the 4th term, 2^3.
%C A331701 The subsequence of odd terms begins: 9, 25, 81, 125, 5041.
%e A331701 9 = 8 + 1.
%e A331701 25 = 16 + 9.
%e A331701 81 = 49 + 32.
%e A331701 125 = 121 + 4.
%e A331701 5041 = 71^2 = 4913 + 128 = 17^3 + 2^7.
%t A331701 Select[#, Last@ # == 1 &][[All, 1]] &@ Fold[Function[{s, k}, Append[s, If[And[! PrimeQ@ k, DivisorSigma[1, k]*EulerPhi[k] > (k - 1)^2], {k, If[AnyTrue[IntegerPartitions[k, {2}], SubsetQ[s[[All, 1]], #] &], 1, 0]}, Nothing]]], {}, Range[10^4]] (* _Michael De Vlieger_, Jan 31 2020 *)
%o A331701 (Python)
%o A331701 from sympy import isprime
%o A331701 TOP = 10**5
%o A331701 primePowers={}
%o A331701 primePowers[1]=1
%o A331701 for x in range(2,TOP):
%o A331701     if isprime(x):
%o A331701         p = pp = x
%o A331701         while pp < TOP**2:
%o A331701             pp *= p
%o A331701             primePowers[pp] = 1
%o A331701 a=[]
%o A331701 pps = sorted(primePowers.keys())[:]
%o A331701 for pp in pps:
%o A331701     for p in pps:
%o A331701         if p*2 > pp: break
%o A331701         if (pp-p) in primePowers:
%o A331701             print(pp)
%o A331701             a.append(pp)
%o A331701             break
%o A331701 print(sorted(a))
%Y A331701 Cf. A025475, A000079.
%K A331701 nonn
%O A331701 1,1
%A A331701 _Alex Ratushnyak_, Jan 25 2020