OFFSET
1,2
COMMENTS
All numbers of the form 3^k2*5^k3*...*p_n^k_n, where k2 >= k3 >= ... >= k_n, sorted.
LINKS
Ray Chandler, Table of n, a(n) for n = 1..10000
David Ryan, Mathematical Harmony Analysis, arXiv preprint arXiv:1603.08904 [cs.SD], 2016-2017.
FORMULA
Sum_{n>=1} 1/a(n) = Product_{n>=2} 1/(1 - 1/A070826(n)) = 1.6241170949... - Amiram Eldar, Oct 20 2020
MATHEMATICA
PrimeExponents[n_] := FactorInteger[n][[All, 2]]; lpe = {}; A147516 = {1}; Do[pe = PrimeExponents[n] // Sort; If[FreeQ[lpe, pe], AppendTo[lpe, pe]; AppendTo[A147516, n]], {n, 3, 40000, 2}]; A147516 (* Jean-François Alcover, Jan 27 2015, after Robert G. Wilson v *)
PROG
(PARI) is(n)=my(k=oo, t); forprime(p=3, , t=valuation(n, p); if(t>k, return(0), k=t); if(k, n/=p^k, return(n==1))) \\ Charles R Greathouse IV, Aug 20 2015
(Python)
from functools import lru_cache
from itertools import count
from sympy import prime, integer_log, primorial
from oeis_sequences.OEISsequences import bisection
def A147516(n):
@lru_cache(maxsize=None)
def g(x, m, j): return sum(g(x//(prime(m)**i), m-1, i) for i in range(j, integer_log(x, prime(m))[0]+1)) if m>2 else max(0, integer_log(x, 3)[0]+1-j)
def f(x):
c = n-1+x
for k in count(2):
if primorial(k)>(x<<1):
break
c -= g(x, k, 1)
return c
return bisection(f, n, n) # Chai Wah Wu, Mar 23 2026
(Python)
from itertools import islice
from heapq import heappop, heappush
from sympy import factorint, prevprime, nextprime
def A147516_gen(): # generator of terms when the first n terms are desired
h, hset = [1], {1}
while True:
yield (m:=heappop(h))
ps = factorint(m)
for p in ps:
if p == 3 or ps[prevprime(p)]>ps[p]:
mp = m*p
if mp not in hset:
heappush(h, mp)
hset.add(mp)
mp = m*nextprime(max(ps.keys(), default=2))
if mp not in hset:
heappush(h, mp)
hset.add(mp)
CROSSREFS
KEYWORD
nonn
AUTHOR
Will Nicholes, Nov 05 2008
EXTENSIONS
Edited and extended by Ray Chandler, Jul 29 2010
STATUS
approved