A069493 6-full numbers: if p divides n then so does p^6.

1, 64, 128, 256, 512, 729, 1024, 2048, 2187, 4096, 6561, 8192, 15625, 16384, 19683, 32768, 46656, 59049, 65536, 78125, 93312, 117649, 131072, 139968, 177147, 186624, 262144, 279936, 373248, 390625, 419904, 524288, 531441, 559872, 746496, 823543, 839808

Offset: 1

Author

Benoit Cloitre, Apr 15 2002

Keywords

• easy
• nonn

Comments

a(m) mod prime(n) > 0 for m < A258603(n); a(A258600(n)) = A030516(n) = prime(n)^6. - Reinhard Zumkeller, Jun 06 2015

Examples

			2^7*3^6 = 93312 is a member (although not of A076470).
		

Links

• Alois P. Heinz, Table of n, a(n) for n = 1..10000

Crossrefs

Cf. A036967, A036966, A001694. Different from A076470.

Cf. A030516, A258603.

Programs

Haskell

import Data.Set (singleton, deleteFindMin, fromList, union)
a069493 n = a069493_list !! (n-1)
a069493_list = 1 : f (singleton z) [1, z] zs where
   f s q6s p6s'@(p6:p6s)
     | m < p6 = m : f (union (fromList $ map (* m) ps) s') q6s p6s'
     | otherwise = f (union (fromList $ map (* p6) q6s) s) (p6:q6s) p6s
     where ps = a027748_row m
           (m, s') = deleteFindMin s
   (z:zs) = a030516_list
-- Reinhard Zumkeller, Jun 03 2015

Mathematica

Join[{1},Select[Range[900000],Min[FactorInteger[#][[All,2]]]>5&]] (* Harvey P. Dale, Mar 03 2018 *)

PARI

for(n=1,560000,if(n*sumdiv(n,d,isprime(d)/d^6)==floor(n*sumdiv(n,d,isprime(d)/d^6)),print1(n,",")))

Python

from math import gcd
from sympy import factorint, integer_nthroot
def A069493(n):
    def f(x):
        c = n+x
        for y1 in range(1,integer_nthroot(x,11)[0]+1):
            if all(d<=1 for d in factorint(y1).values()):
                for y2 in range(1,integer_nthroot(z2:=x//y1**11,10)[0]+1):
                    if gcd(y2,y1)==1 and all(d<=1 for d in factorint(y2).values()):
                        for y3 in range(1,integer_nthroot(z3:=z2//y2**10,9)[0]+1):
                            if gcd(y3,y1)==1 and gcd(y3,y2)==1 and all(d<=1 for d in factorint(y3).values()):
                                for y4 in range(1,integer_nthroot(z4:=z3//y3**9,8)[0]+1):
                                    if gcd(y4,y1)==1 and gcd(y4,y2)==1 and gcd(y4,y3)==1 and all(d<=1 for d in factorint(y4).values()):
                                        for y5 in range(1,integer_nthroot(z5:=z4//y4**8,7)[0]+1):
                                            if gcd(y5,y1)==1 and gcd(y5,y2)==1 and gcd(y5,y3)==1 and gcd(y5,y4)==1 and all(d<=1 for d in factorint(y5).values()):
                                                c -= integer_nthroot(z5//y5**7,6)[0]
        return c
    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
    return bisection(f,n,n) # Chai Wah Wu, Sep 10 2024

Formula

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