A303661 Powers of squarefree semiprimes that are not squarefree.

36, 100, 196, 216, 225, 441, 484, 676, 1000, 1089, 1156, 1225, 1296, 1444, 1521, 2116, 2601, 2744, 3025, 3249, 3364, 3375, 3844, 4225, 4761, 5476, 5929, 6724, 7225, 7396, 7569, 7776, 8281, 8649, 8836, 9025, 9261, 10000, 10648, 11236, 12321, 13225, 13924, 14161, 14884

Offset: 1

Author

Ilya Gutkovskiy, Apr 28 2018

Keywords

• nonn

Examples

			1089 is in the sequence because 1089 = 3^2*11^2.
1296 is in the sequence because 1296 = 2^4*3^4.
		

Links

• Amiram Eldar, Table of n, a(n) for n = 1..10000
• Eric Weisstein's World of Mathematics, Squarefree.
• Eric Weisstein's World of Mathematics, Semiprime.

Crossrefs

Cf. A006881, A013929, A072774, A072777, A085155, A123711, A200511, A246547, A303606.

Programs

Mathematica

Select[Range[15000], Length[Union[FactorInteger[#][[All, 2]]]] == 1 && PrimeNu[#] == 2 && ! SquareFreeQ[#] &]
seq[max_] := Module[{sp = Select[Range[Floor@Sqrt[max]], SquareFreeQ[#] && PrimeNu[#] == 2 &], s = {}}, Do[s = Join[s, sp[[k]]^Range[2, Floor@Log[sp[[k]], max]]], {k, 1, Length[sp]}]; Union@s]; seq[10000] (* Amiram Eldar, Feb 12 2021 *)

Python

from math import isqrt
from sympy import primepi, primerange, integer_nthroot
def A303661(n):
    def g(x): return int(-(t:=primepi(s:=isqrt(x)))-(t*(t-1)>>1)+sum(primepi(x//k) for k in primerange(1, s+1)))
    def f(x): return n-1+x-sum(g(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length()))
    kmin, kmax = 1,2
    while f(kmax) >= kmax:
        kmax <<= 1
    while True:
        kmid = kmax+kmin>>1
        if f(kmid) < kmid:
            kmax = kmid
        else:
            kmin = kmid
        if kmax-kmin <= 1:
            break
    return kmax # Chai Wah Wu, Aug 19 2024

Formula

Sum_{n>=1} 1/a(n) = Sum_{n>=1} 1/((A006881(n)-1)*A006881(n)) = Sum_{k>=2} (P(k)^2 - P(2*k))/2 = 0.07160601536406295068..., where P(k) is the prime zeta function. - Amiram Eldar, Feb 12 2021