A069492 5-full numbers: if a prime p divides k then so does p^5.

1, 32, 64, 128, 243, 256, 512, 729, 1024, 2048, 2187, 3125, 4096, 6561, 7776, 8192, 15552, 15625, 16384, 16807, 19683, 23328, 31104, 32768, 46656, 59049, 62208, 65536, 69984, 78125, 93312, 100000, 117649, 124416, 131072, 139968, 161051

Offset: 1

Author

Benoit Cloitre, Apr 15 2002

Keywords

• easy
• nonn

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)

Crossrefs

Cf. A036967, A036966, A001694.

Cf. A050997.

Cf. A258602.

Programs

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

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

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

Formula

Sum_{n>=1} 1/a(n) = Product_{p prime} (1 + 1/(p^4*(p-1))) = 1.0695724994489739263413712783666538355049945684326048537289707764272637... - Amiram Eldar, Jul 09 2020