OFFSET
1,2
COMMENTS
If p|n then at least p^4|n.
Subsequence of A036967. - R. J. Mathar, May 27 2011
LINKS
Reinhard Zumkeller, Table of n, a(n) for n = 1..10000
FORMULA
Sum_{n>=1} 1/a(n) = 3 - zeta(2) - zeta(3) + Sum_{k>=2} mu(k)*(3 - zeta(k) - zeta(2*k) - zeta(3*k)) = 1.1473274274... . - Amiram Eldar, Dec 03 2022
MATHEMATICA
a = {1}; Do[ If[ Apply[ GCD, Last[ Transpose[ FactorInteger[n]]]] > 3, a = Append[a, n]; Print[n]], {n, 2, 131071}]; a
PROG
(Haskell)
import qualified Data.Set as Set (null)
import Data.Set (empty, insert, deleteFindMin)
a076468 n = a076468_list !! (n-1)
a076468_list = 1 : f [2..] empty where
f xs'@(x:xs) s | Set.null s || m > x ^ 4 = f xs $ insert (x ^ 4, x) s
| m == x ^ 4 = f xs s
| otherwise = m : f xs' (insert (m * b, b) s')
where ((m, b), s') = deleteFindMin s
-- Reinhard Zumkeller, Jun 19 2013
(Python)
from sympy import mobius, integer_nthroot
def A076468(n):
def f(x): return int(n+2+x-integer_nthroot(x, 4)[0]-(integer_nthroot(x, 6)[0]<<1)-integer_nthroot(x, 9)[0]+sum(mobius(k)*(integer_nthroot(x, k)[0]+integer_nthroot(x, k<<1)[0]+integer_nthroot(x, 3*k)[0]-3) for k in range(5, 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 14 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
Robert G. Wilson v, Oct 14 2002
STATUS
approved