cp's OEIS Frontend

A036967 4-full numbers: if a prime p divides k then so does p^4.

%I A036967 #49 Sep 11 2024 00:31:56
%S A036967 1,16,32,64,81,128,243,256,512,625,729,1024,1296,2048,2187,2401,2592,
%T A036967 3125,3888,4096,5184,6561,7776,8192,10000,10368,11664,14641,15552,
%U A036967 15625,16384,16807,19683,20000,20736,23328,28561,31104,32768,34992
%N A036967 4-full numbers: if a prime p divides k then so does p^4.
%C A036967 a(m) mod prime(n) > 0 for m < A258601(n); a(A258601(n)) = A030514(n) = prime(n)^4. - _Reinhard Zumkeller_, Jun 06 2015
%D A036967 E. Kraetzel, Lattice Points, Kluwer, Chap. 7, p. 276.
%H A036967 Alois P. Heinz, <a href="/A036967/b036967.txt">Table of n, a(n) for n = 1..10000</a> (first 300 terms from T. D. Noe)
%F A036967 Sum_{n>=1} 1/a(n) = Product_{p prime} (1 + 1/(p^3*(p-1))) = 1.1488462139214317030108176090790939019972506733993367867997411290952527... - _Amiram Eldar_, Jul 09 2020
%t A036967 Join[{1},Select[Range[35000],Min[Transpose[FactorInteger[#]][[2]]]>3&]] (* _Harvey P. Dale_, Jun 05 2012 *)
%o A036967 (Haskell)
%o A036967 import Data.Set (singleton, deleteFindMin, fromList, union)
%o A036967 a036967 n = a036967_list !! (n-1)
%o A036967 a036967_list = 1 : f (singleton z) [1, z] zs where
%o A036967    f s q4s p4s'@(p4:p4s)
%o A036967      | m < p4 = m : f (union (fromList $ map (* m) ps) s') q4s p4s'
%o A036967      | otherwise = f (union (fromList $ map (* p4) q4s) s) (p4:q4s) p4s
%o A036967      where ps = a027748_row m
%o A036967            (m, s') = deleteFindMin s
%o A036967    (z:zs) = a030514_list
%o A036967 -- _Reinhard Zumkeller_, Jun 03 2015
%o A036967 (PARI) is(n)=n==1 || vecmin(factor(n)[,2])>3 \\ _Charles R Greathouse IV_, Sep 17 2015
%o A036967 (PARI)
%o A036967 M(v,u,lim)=vecsort(concat(vector(#v, i, my(m=lim\v[i]); v[i]*select(t->t<=m, u))))
%o A036967 Gen(lim,k)={my(v=[1]); forprime(p=2, sqrtnint(lim, k), v=M(v, concat([1], vector(logint(lim,p)-k+1,e,p^(e+k-1))), lim));v}
%o A036967 Gen(35000,4) \\ _Andrew Howroyd_, Sep 10 2024
%o A036967 (Python)
%o A036967 from sympy import factorint
%o A036967 A036967_list = [n for n in range(1,10**5) if min(factorint(n).values(),default=4) >= 4] # _Chai Wah Wu_, Aug 18 2021
%o A036967 (Python)
%o A036967 from math import gcd
%o A036967 from sympy import integer_nthroot, factorint
%o A036967 def A036967(n):
%o A036967     def bisection(f,kmin=0,kmax=1):
%o A036967         while f(kmax) > kmax: kmax <<= 1
%o A036967         while kmax-kmin > 1:
%o A036967             kmid = kmax+kmin>>1
%o A036967             if f(kmid) <= kmid:
%o A036967                 kmax = kmid
%o A036967             else:
%o A036967                 kmin = kmid
%o A036967         return kmax
%o A036967     def f(x):
%o A036967         c = n+x
%o A036967         for u in range(1,integer_nthroot(x,7)[0]+1):
%o A036967             if all(d<=1 for d in factorint(u).values()):
%o A036967                 for w in range(1,integer_nthroot(a:=x//u**7,6)[0]+1):
%o A036967                     if gcd(w,u)==1 and all(d<=1 for d in factorint(w).values()):
%o A036967                         for y in range(1,integer_nthroot(z:=a//w**6,5)[0]+1):
%o A036967                             if gcd(w,y)==1 and gcd(u,y)==1 and all(d<=1 for d in factorint(y).values()):
%o A036967                                 c -= integer_nthroot(z//y**5,4)[0]
%o A036967         return c
%o A036967     return bisection(f,n,n) # _Chai Wah Wu_, Sep 10 2024
%Y A036967 A030514 is a subsequence.
%Y A036967 Cf. A001694, A036966, A046101, A258601.
%K A036967 easy,nonn,nice
%O A036967 1,2
%A A036967 _N. J. A. Sloane_
%E A036967 More terms from _Erich Friedman_
%E A036967 Corrected by _Vladeta Jovovic_, Aug 17 2002