A280609 - OEIS

%I #14 Sep 13 2024 06:54:31

%S 9,25,27,49,121,125,169,243,289,343,361,529,841,961,1331,1369,1681,

%T 1849,2187,2197,2209,2809,3125,3481,3721,4489,4913,5041,5329,6241,

%U 6859,6889,7921,9409,10201,10609,11449,11881,12167,12769,16129,16807,17161,18769,19321,22201,22801,24389,24649,26569,27889,29791,29929

%N Odd prime powers with prime exponents.

%C Intersection of A053810 and A061345.

%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/PrimePower.html">Prime Power</a>.

%F a(n) = p^q, where p, q are primes and p > 2.

%F Sum_{n>=1} 1/a(n) = Sum_{p prime} P(p) - A051006 = 0.25699271237062131298..., where P(s) is the prime zeta function. - _Amiram Eldar_, Sep 13 2024

%e 9 is in the sequence because 9 = 3^2;

%e 25 is in the sequence because 25 = 5^2;

%e 27 is in the sequence because 27 = 3^3, etc.

%t Select[Range[30000], PrimePowerQ[#1] && PrimeQ[PrimeOmega[#1]] && Mod[#1, 2] == 1 & ]

%o (Python)

%o from sympy import primepi, integer_nthroot, primerange

%o def A280609(n):

%o     def f(x): return int(n+x-sum(primepi(integer_nthroot(x, p)[0])-1 for p in primerange(x.bit_length())))

%o     def bisection(f,kmin=0,kmax=1):

%o         while f(kmax) > kmax: kmax <<= 1

%o         while kmax-kmin > 1:

%o             kmid = kmax+kmin>>1

%o             if f(kmid) <= kmid:

%o                 kmax = kmid

%o             else:

%o                 kmin = kmid

%o         return kmax

%o     return bisection(f,n,n) # _Chai Wah Wu_, Sep 12 2024

%Y Cf. A000961, A034785, A051006, A053810, A061345, A078422, A246547, A246551, A246655.

%K nonn,easy

%O 1,1

%A _Ilya Gutkovskiy_, Jan 06 2017