OFFSET
1,1
LINKS
Michael S. Branicky, Table of n, a(n) for n = 1..10001 (first 1001 terms from Harvey P. Dale)
FORMULA
a(n) = 3*A002473(n). - Chai Wah Wu, Sep 18 2024
Sum_{n>=1} 1/a(n) = 35/24. - Amiram Eldar, Sep 23 2024
MATHEMATICA
Select[3*Range[200], FactorInteger[#][[-1, 1]]<10&] (* Harvey P. Dale, Apr 10 2019 *)
PROG
(Python)
from sympy import integer_log
def A085126(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 i in range(integer_log(x, 7)[0]+1):
for j in range(integer_log(m:=x//7**i, 5)[0]+1):
for k in range(integer_log(r:=m//5**j, 3)[0]+1):
c -= (r//3**k).bit_length()
return c
return bisection(f, n, n)*3 # Chai Wah Wu, Sep 17 2024
(Python) # faster for initial segment of sequence
import heapq
from itertools import islice
def A085126gen(): # generator of terms
v, oldv, h, psmooth_primes, = 1, 0, [1], [2, 3, 5, 7]
while True:
v = heapq.heappop(h)
if v != oldv:
yield 3*v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
print(list(islice(A085126gen(), 65))) # Michael S. Branicky, Sep 17 2024
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
Amarnath Murthy, Jul 06 2003
EXTENSIONS
More terms from David Wasserman, Jan 28 2005
Offset changed to 1 by Michael S. Branicky, Sep 17 2024
STATUS
approved