OFFSET
1,2
COMMENTS
Let i = a(n-2), j = a(n-1). For k > 1, m >= 1, a(n) = m*prime(k) iff rad(i*j) = primorial(k-1), and this is the m-th such occurrence. This suggests the late appearance of most primes (namely those >= 7), apparent in the lowest part of scatterplot, where for example a(717126), a(63056215) = 31, 37 respectively.
As Scott R. Shannon has just observed, the following proof is incomplete, since it requires a proof that every even number appears. Even the induction step seems a little dubious. - N. J. A. Sloane, Mar 18 2023
All multiples of all primes appear in the sequence, for if not there is a least prime p such that m*p is not a term for any [some?] m >= 1. Choose any prime q < p; then every multiple of q must appear, so then p*q must be a term; contradiction since this is a multiple of p. [But what if p = 2?]
Corollary: This sequence is a permutation of the positive integers. [This question appears to be still open. - N. J. A. Sloane, Mar 18 2023]
Conjecture: The primes appear in their natural order.
LINKS
Winston de Greef, Table of n, a(n) for n = 1..10000
Michael De Vlieger, Log log scatterplot of a million terms showing primes in red.
Michael De Vlieger, Log log scatterplot of a(n), n = 1..2^16, showing primes in red, composite prime powers in gold, squarefree composites in green, and numbers neither prime power nor squarefree in blue.
EXAMPLE
a(3) must be 3 because a(1,2) = 1,2 and 3 is the least prime which does not divide 2.
a(4) = 5 since this is the least multiple of the smallest prime which does not divide 2*3 = 6.
a(8) = 7 because a(6,7) = 6,10 and 7 is the smallest prime which does not divide 60, rad(60) = 2*3*5 = 30.
a(19,20) = 18,28, and 5 is the smallest prime not dividing rad(18*28) = 42. Since multiples of 5 have appeared 5 times already, a(20) = 6*5 = 30.
MAPLE
R:= 1, 2: S:= {1, 2}:
for i from 3 to 100 do
s:= R[i-2]*R[i-1]:
p:= 2;
while s mod p = 0 do p:= nextprime(p) od:
for r from p by p while member(r, S) do od:
R:= R, r; S:= S union {r}
od:
R; # Robert Israel, Mar 08 2023
MATHEMATICA
nn = 2^10; c[_] = False; q[_] = 1;
Array[Set[{a[#], c[#]}, {#, True}] &, 2];
Set[{i, j}, {a[1], a[2]}]; u = 3;
Do[(k = q[#];
While[c[k #], k++]; k *= #;
While[c[# q[#]], q[#]++]) &[(p = 2;
While[Divisible[i j, p], p = NextPrime[p]]; p)];
Set[{a[n], c[k], i, j}, {k, True, j, k}];
If[k == u, While[c[u], u++]], {n, 3, nn}];
Array[a, nn] (* Michael De Vlieger, Mar 08 2023 *)
PROG
(PARI) findp(n) = forprime(p=2, , if (n%p, return(p)));
lista(nn) = my(va = vector(nn, k, if (k<=2, k))); for (n=3, nn, my(vsa = vecsort(va), p=findp(va[n-1]*va[n-2]), k=p); while (vecsearch(vsa, k), k+=p); va[n] = k; ); va; \\ Michel Marcus, Mar 09 2023
(Python)
from itertools import count, islice
from sympy import prime, primefactors, primepi
def A359804_gen(): # generator of terms
aset, bset, cset = set(), {1}, {1, 2}
yield from (1, 2)
while True:
for i in count(1):
if not (i in aset or i in bset):
p = prime(i)
for j in count(1):
if (m:=j*p) not in cset:
yield m
cset.add(m)
break
break
aset, bset = bset, set(map(primepi, primefactors(m)))
CROSSREFS
KEYWORD
nonn
AUTHOR
David James Sycamore, Mar 08 2023
STATUS
approved