OFFSET
1,1
COMMENTS
This sequence is the first step toward candidates for odd perfect numbers, A228058.
LINKS
T. D. Noe, Table of n, a(n) for n = 1..10000
Raghavendra Bhat, Distribution of Square Prime Numbers, arXiv:2109.10238 [math.NT], 2021.
Raghavendra Bhat, An Algebraic Structure for Square-Prime Numbers, arXiv:2303.14296 [math.GM], 2023.
Raghavendra Bhat, Cristian Cobeli, and Alexandru Zaharescu, Filtered rays over iterated absolute differences on layers of integers, arXiv:2309.03922 [math.NT], 2023. See p. 9 (Sect. 3.1).
David Harvey, Faster enumeration of primes, arXiv:2606.22851 [math.NT], 2026. See p. 12 (Remark 1.7).
FORMULA
Bhat proves there are ~ (Pi^2/6-1)*x/log x members of this sequence up to x, so a(n) ~ kn log n with k = 6/(Pi^2-6) = 1.550546.... - Charles R Greathouse IV, Oct 01 2021
MATHEMATICA
nn = 300; Union[Select[Flatten[Table[p*n^2, {p, Prime[Range[PrimePi[nn/4]]]}, {n, 2, Sqrt[nn/2]}]], # < nn &]]
PROG
(Haskell)
import Data.List (partition)
a228056 n = a228056_list !! (n-1)
a228056_list = filter f [1..] where
f x = length us == 1 && (head us > 1 || not (null vs)) where
(us, vs) = partition odd $ a124010_row x
-- Reinhard Zumkeller, Aug 14 2013
(PARI) list(lim)=my(v=List()); forfactored(n=2, lim\1, my(e=n[2][, 2]); if(vecsum(e%2)==1 && e!=[1]~, listput(v, n[1]))); Vec(v); \\ Charles R Greathouse IV, Oct 01 2021
(Python)
from math import isqrt
from sympy import primepi
def A228056(n):
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
kmin = kmax >> 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x): return n+x-sum(primepi(x//y**2) for y in range(2, isqrt(x)+1))
return bisection(f, n, n) # Chai Wah Wu, Jun 06 2025
CROSSREFS
KEYWORD
nonn
AUTHOR
T. D. Noe, Aug 13 2013
STATUS
approved