A359804 a(1) = 1, a(2) = 2; thereafter let p be the smallest prime that does not divide a(n-2)*a(n-1), then a(n) is the smallest multiple of p that is not yet in the sequence.

1, 2, 3, 5, 4, 6, 10, 7, 9, 8, 15, 14, 11, 12, 20, 21, 22, 25, 18, 28, 30, 33, 35, 16, 24, 40, 42, 44, 45, 49, 26, 27, 50, 56, 36, 55, 63, 32, 60, 70, 66, 13, 65, 34, 39, 75, 38, 77, 48, 80, 84, 88, 85, 51, 46, 90, 91, 99, 52, 95, 54, 98, 100, 57, 105, 58, 110, 69, 112, 115, 72, 119, 120, 121

Offset: 1

Author

David James Sycamore, Mar 08 2023

Keywords

• nonn

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.

Examples

			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.
		

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.

Crossrefs

Cf. A002110, A007947, A053669.

See also A361502, A361503, A361504, A361505.

A351495 has a very similar definition.

Programs

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 *)

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)))
A359804_list = list(islice(A359804_gen(),30)) # Chai Wah Wu, Mar 18 2023