OFFSET
1,2
COMMENTS
a(m) mod prime(n) > 0 for m < A258602(n); a(A258602(n)) = A050997(n) = prime(n)^5. - Reinhard Zumkeller, Jun 06 2015
LINKS
Reinhard Zumkeller, Table of n, a(n) for n = 1..10000 (terms n=148..10000 corrected by Andrew Howroyd)
FORMULA
Sum_{n>=1} 1/a(n) = Product_{p prime} (1 + 1/(p^4*(p-1))) = 1.0695724994489739263413712783666538355049945684326048537289707764272637... - Amiram Eldar, Jul 09 2020
PROG
(PARI) for(n=1, 250000, if(n*sumdiv(n, d, isprime(d)/d^5)==floor(n*sumdiv(n, d, isprime(d)/d^5)), print1(n, ", ")))
(PARI) \\ Gen(limit, k) defined in A036967.
Gen(170000, 5) \\ Andrew Howroyd, Sep 10 2024
(Haskell)
import Data.Set (singleton, deleteFindMin, fromList, union)
a069492 n = a069492_list !! (n-1)
a069492_list = 1 : f (singleton z) [1, z] zs where
f s q5s p5s'@(p5:p5s)
| m < p5 = m : f (union (fromList $ map (* m) ps) s') q5s p5s'
| otherwise = f (union (fromList $ map (* p5) q5s) s) (p5:q5s) p5s
where ps = a027748_row m
(m, s') = deleteFindMin s
(z:zs) = a050997_list
-- Reinhard Zumkeller, Jun 03 2015
(Python)
from math import gcd
from sympy import integer_nthroot, factorint
def A069492(n):
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c = n+x
for t in range(1, integer_nthroot(x, 9)[0]+1):
if all(d<=1 for d in factorint(t).values()):
for u in range(1, integer_nthroot(s:=x//t**9, 8)[0]+1):
if gcd(t, u)==1 and all(d<=1 for d in factorint(u).values()):
for w in range(1, integer_nthroot(a:=s//u**8, 7)[0]+1):
if gcd(u, w)==1 and gcd(t, w)==1 and all(d<=1 for d in factorint(w).values()):
for y in range(1, integer_nthroot(z:=a//w**7, 6)[0]+1):
if gcd(w, y)==1 and gcd(u, y)==1 and gcd(t, y)==1 and all(d<=1 for d in factorint(y).values()):
c -= integer_nthroot(z//y**6, 5)[0]
return c
return bisection(f, n, n) # Chai Wah Wu, Sep 10 2024
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
Benoit Cloitre, Apr 15 2002
STATUS
approved