OFFSET
1,2
COMMENTS
Also, distinct terms appearing in A352598, or terms of the form 4^i * 9^j * 25^k * 49^m for i, j, k, m >= 0.
LINKS
Michael De Vlieger, Table of n, a(n) for n = 1..10000
FORMULA
a(n) = A002473(n)^2. - Pontus von Brömssen, Mar 24 2022
Sum_{n>=1} 1/a(n) = 1225/768. - Amiram Eldar, Mar 24 2022
EXAMPLE
49 = 7*7, 81 = (3*3)*(3*3), and 100 = (2*5)*(2*5) are terms.
MATHEMATICA
Select[Range[120], Max[FactorInteger[#][[;; , 1]]] <= 7 &]^2 (* Amiram Eldar, Mar 24 2022 *)
With[{n = 15000}, Union@ Flatten@ Table[2^(2 a)*3^(2 b)*5^(2 c)*7^(2 d), {a, 0, Log[4, n]}, {b, 0, Log[9, n/(2^(2 a))]}, {c, 0, Log[25, n/(2^(2 a)*3^(2 b))]}, {d, 0, Log[49, n/(2^(2 a)*3^(2 b)*5^(2 c))]}]] (* Michael De Vlieger, Mar 26 2022 *)
PROG
(Python)
from itertools import count, islice
def agen():
for i in count(1):
k = i
for p in [2, 3, 5, 7]:
while k%p == 0:
k //= p
if k == 1:
yield i*i
print(list(islice(agen(), 50)))
(Python)
from sympy import integer_log
def A352618(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)**2 # Chai Wah Wu, Sep 17 2024
(Python) # faster for initial segment of sequence
import heapq
from itertools import islice
from sympy import primerange
def A352618gen(p=7): # generator of 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*v
oldv = v
for p in psmooth_primes:
heapq.heappush(h, v*p)
print(list(islice(A352618gen(), 65))) # Michael S. Branicky, Sep 17 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
Michael S. Branicky, Mar 24 2022
STATUS
approved