OFFSET
1,2
COMMENTS
a(m) mod prime(n) > 0 for m < A258601(n); a(A258601(n)) = A030514(n) = prime(n)^4. - Reinhard Zumkeller, Jun 06 2015
REFERENCES
E. Kraetzel, Lattice Points, Kluwer, Chap. 7, p. 276.
LINKS
Alois P. Heinz, Table of n, a(n) for n = 1..10000 (first 300 terms from T. D. Noe)
FORMULA
Sum_{n>=1} 1/a(n) = Product_{p prime} (1 + 1/(p^3*(p-1))) = 1.1488462139214317030108176090790939019972506733993367867997411290952527... - Amiram Eldar, Jul 09 2020
MATHEMATICA
Join[{1}, Select[Range[35000], Min[Transpose[FactorInteger[#]][[2]]]>3&]] (* Harvey P. Dale, Jun 05 2012 *)
PROG
(Haskell)
import Data.Set (singleton, deleteFindMin, fromList, union)
a036967 n = a036967_list !! (n-1)
a036967_list = 1 : f (singleton z) [1, z] zs where
f s q4s p4s'@(p4:p4s)
| m < p4 = m : f (union (fromList $ map (* m) ps) s') q4s p4s'
| otherwise = f (union (fromList $ map (* p4) q4s) s) (p4:q4s) p4s
where ps = a027748_row m
(m, s') = deleteFindMin s
(z:zs) = a030514_list
-- Reinhard Zumkeller, Jun 03 2015
(PARI) is(n)=n==1 || vecmin(factor(n)[, 2])>3 \\ Charles R Greathouse IV, Sep 17 2015
(PARI)
M(v, u, lim)=vecsort(concat(vector(#v, i, my(m=lim\v[i]); v[i]*select(t->t<=m, u))))
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}
Gen(35000, 4) \\ Andrew Howroyd, Sep 10 2024
(Python)
from sympy import factorint
A036967_list = [n for n in range(1, 10**5) if min(factorint(n).values(), default=4) >= 4] # Chai Wah Wu, Aug 18 2021
(Python)
from math import gcd
from sympy import integer_nthroot, factorint
def A036967(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 u in range(1, integer_nthroot(x, 7)[0]+1):
if all(d<=1 for d in factorint(u).values()):
for w in range(1, integer_nthroot(a:=x//u**7, 6)[0]+1):
if gcd(w, u)==1 and all(d<=1 for d in factorint(w).values()):
for y in range(1, integer_nthroot(z:=a//w**6, 5)[0]+1):
if gcd(w, y)==1 and gcd(u, y)==1 and all(d<=1 for d in factorint(y).values()):
c -= integer_nthroot(z//y**5, 4)[0]
return c
return bisection(f, n, n) # Chai Wah Wu, Sep 10 2024
CROSSREFS
KEYWORD
easy,nonn,nice
AUTHOR
EXTENSIONS
More terms from Erich Friedman
Corrected by Vladeta Jovovic, Aug 17 2002
STATUS
approved