OFFSET
1,1
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000
Eric Weisstein's World of Mathematics, Prime Power.
Eric Weisstein's World of Mathematics, Squarefree.
FORMULA
Sum_{n>=1} 1/a(n) = Sum_{n>=1} 1/((A120944(n)-1)*A120944(n)) = Sum_{k>=2} (zeta(k)/zeta(2*k) - P(k) - 1) = 0.07547719891508850482..., where P(k) is the prime zeta function. - Amiram Eldar, Feb 12 2021
EXAMPLE
196 is in the sequence because 196 = 2^2*7^2.
4900 is in the sequence because 4900 = 2^2*5^2*7^2.
MATHEMATICA
Select[Range[12000], Length[Union[FactorInteger[#][[All, 2]]]] == 1 && ! SquareFreeQ[#] && ! PrimePowerQ[#] &]
seq[max_] := Module[{sp = Select[Range[Floor@Sqrt[max]], SquareFreeQ[#] && PrimeNu[#] > 1 &], s = {}}, Do[s = Join[s, sp[[k]]^Range[2, Floor@Log[sp[[k]], max]]], {k, 1, Length[sp]}]; Union@s]; seq[10^4] (* Amiram Eldar, Feb 12 2021 *)
PROG
(Python)
from math import isqrt
from sympy import mobius, primepi, integer_nthroot
def A303606(n):
def g(x): return int(sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))-primepi(x))
def f(x): return n-3+x+(y:=x.bit_length())-sum(g(integer_nthroot(x, k)[0]) for k in range(2, y))
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
CROSSREFS
KEYWORD
nonn
AUTHOR
Ilya Gutkovskiy, Apr 26 2018
STATUS
approved