cp's OEIS Frontend

A303661 Powers of squarefree semiprimes that are not squarefree.

%I A303661 #15 Feb 16 2025 08:33:54
%S A303661 36,100,196,216,225,441,484,676,1000,1089,1156,1225,1296,1444,1521,
%T A303661 2116,2601,2744,3025,3249,3364,3375,3844,4225,4761,5476,5929,6724,
%U A303661 7225,7396,7569,7776,8281,8649,8836,9025,9261,10000,10648,11236,12321,13225,13924,14161,14884
%N A303661 Powers of squarefree semiprimes that are not squarefree.
%H A303661 Amiram Eldar, <a href="/A303661/b303661.txt">Table of n, a(n) for n = 1..10000</a>
%H A303661 Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/Squarefree.html">Squarefree</a>.
%H A303661 Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/Semiprime.html">Semiprime</a>.
%F A303661 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
%e A303661 1089 is in the sequence because 1089 = 3^2*11^2.
%e A303661 1296 is in the sequence because 1296 = 2^4*3^4.
%t A303661 Select[Range[15000], Length[Union[FactorInteger[#][[All, 2]]]] == 1 && PrimeNu[#] == 2 && ! SquareFreeQ[#] &]
%t A303661 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 *)
%o A303661 (Python)
%o A303661 from math import isqrt
%o A303661 from sympy import primepi, primerange, integer_nthroot
%o A303661 def A303661(n):
%o A303661     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)))
%o A303661     def f(x): return n-1+x-sum(g(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length()))
%o A303661     kmin, kmax = 1,2
%o A303661     while f(kmax) >= kmax:
%o A303661         kmax <<= 1
%o A303661     while True:
%o A303661         kmid = kmax+kmin>>1
%o A303661         if f(kmid) < kmid:
%o A303661             kmax = kmid
%o A303661         else:
%o A303661             kmin = kmid
%o A303661         if kmax-kmin <= 1:
%o A303661             break
%o A303661     return kmax # _Chai Wah Wu_, Aug 19 2024
%Y A303661 Cf. A006881, A013929, A072774, A072777, A085155, A123711, A200511, A246547, A303606.
%K A303661 nonn
%O A303661 1,1
%A A303661 _Ilya Gutkovskiy_, Apr 28 2018