Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).
reviewed
approved
proposed
reviewed
editing
proposed
return bisection(f, n, n) # Chai Wah Wu, Sep 16 2024
(Python)
from sympy import integer_log, prevprime
def A051038(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 g(x, m): return sum((x//3**i).bit_length() for i in range(integer_log(x, 3)[0]+1)) if m==3 else sum(g(x//(m**i), prevprime(m))for i in range(integer_log(x, m)[0]+1))
def f(x): return n+x-g(x, 11)
return bisection(f, n, n) # Chai Wah Wu, Sep 16 2024
approved
editing
reviewed
approved
proposed
reviewed
editing
proposed
print(list(islice(agen(), 6567))) # Michael S. Branicky, Nov 20 2022
(Python)
import heapq
from itertools import islice
from sympy import primerange
def agen(p=11): # generate all p-smooth terms
v, oldv, h, psmooth_primes, = 1, 0, [1], list(primerange(1, p+1))
while True:
v = heapq.heappop(h)
if v != oldv:
yield v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
print(list(islice(agen(), 65))) # Michael S. Branicky, Nov 20 2022
approved
editing