A358093 - OEIS

%I #65 Nov 14 2022 00:34:30

%S 1,2,3,5,4,7,9,10,13,23,6,29,35,8,37,39,38,77,115,12,11,17,14,31,15,

%T 46,61,107,42,149,191,170,19,21,20,961,41,18,47,53,40,63,29791,26,57,

%U 83,70,51,121,62,73,45,22,1369,59,24,65,71,34,105,139,122,87,209,74

%N a(n) = n for 1 <= n <= 2; thereafter a(n) is the least unused m such that rad(m) = rad(rad(a(n-1)) + rad(a(n-2))), where rad(m) = A007947(m).

%C In other words, a(n) is the least m having the same squarefree kernel as the sum of the squarefree kernels of a(n-1) and a(n-2).

%C Primes do not appear in natural order.

%C Conjectured to be a permutation of the positive integers.

%C From _Michael De Vlieger_, Nov 09 2022: (Start)

%C Let i = rad(a(n-2)), j = rad(a(n-1)), and s = rad(i+j). If s = p (prime), then a(n) = m = p^e, e >= 1. Because a(n) = m such that rad(m) = p, the j-th occasion of s = p implies a(n) = p^j.

%C By the same token, generally, the j-th occasion of s implies a(n) = M*s, where M is regular to s, i.e., M is a product restricted to primes p|s, with M appearing in order. For example, if s = 10, then the j-th occasion of s = 10 implies a(n) = 10*A003592(j). These are consequences of conservation of squarefree kernel s, the lexical and greedy nature of the sequence.

%C It is clear that a(n) >= s, since rad(a(n)) = s, hence s | a(n), and on account of the lexical and greedy nature of the sequence.

%C Prime a(n) implies s = a(n), while s < a(n) for a(n) that are composite prime powers (i.e., a(n) in A246547) since s is in those cases prime.(End)

%C From _Michael De Vlieger_, Nov 13 2022: (Start)

%C Adjacent terms i and j are coprime as consequence of s mod p = i mod p + j mod p and a(n) = M*s such that rad(M) = rad(s). p | a(n-2) and p | a(n-1) implies p | a(n), and subsequent adjacent terms are likewise divisible by p, contradicting initial a(1) coprime to a(2). Because the sequence starts with dissimilar residues (mod p), p divides only one of {a(n-2), a(n-1), a(n)}.

%C 2 | a(n) for n mod 3 = 2. (End)

%H Michael De Vlieger, <a href="/A358093/b358093.txt">Table of n, a(n) for n = 1..1500</a>

%H Michael De Vlieger, <a href="/A358093/a358093.png">Log log scatterplot of a(n)</a>, n = 1..1500, highlighting primes in red, other prime powers in gold, and other numbers in blue.

%e a(5) = 4 because rad(rad(3) + rad(5)) = rad(3 + 5) = rad(8) = 2, and 4 is the least unused number m such that rad(m) = 2.

%e a(36) = rad(rad(a(34)) + rad(a(35)) = rad(rad(21) + rad(20)) = rad(31) = 31, and since 31 has appeared at a(24), a(36) = 31^2 = 961.

%e a(43) = 29791 (31^3) because rad(40) + rad(63) = 31; we already have 31 and 31^2.

%t nn = 65; c[_] = False; p[_] = q[_] = 1; Array[Set[{a[#], c[#]}, {#, True}] &, 3]; Array[(q[#]++; p[#]++) &[Times @@ FactorInteger[a[#]][[All, 1]]] &, 2, 2]; u = 4; Set[{i, j}, Map[Times @@ FactorInteger[#][[All, 1]] &, {a[2], a[3]}]]; Do[k = Times @@ FactorInteger[i + j][[All, 1]]; If[PrimeQ[k], m = k^p[k]; p[k]++, m = q[k]; While[Nand[! c[m k], PowerMod[k, k, m] == 0], m++]; m *= k]; q[k]++; Set[{a[n], c[m], i, j}, {m, True, j, k}]; If[m == u, While[c[u], u++]], {n, 4, nn}]; Array[a, nn] (* _Michael De Vlieger_, Nov 09 2022 *)

%o (PARI) first(n) = { n = max(n, 3); res = vector(n); for(i = 1, 3, res[i] = i); for(i = 4, n, target = rad(rad(res[i-1]) + rad(res[i-2])); resset = Set(res); pr = factor(target)[, 1]; step = target; if(omega(target) >= 1, forstep(j = target, oo, step, if(rad(j) == target, if(#setminus(Set(j), resset) == 1, res[i] = j; next(2) ) ) ) , for(j = 1, oo, if(#setminus(Set(target ^ j), resset) == 1, res[i] = j; next(2) ) ) ) ); res }

%o rad(n) = factorback(factorint(n)[, 1]); \\ _David A. Corneth_, Nov 09 2022

%o (PARI) nn=300; a=List([1,2]); rad(n)= factorback(factorint(n)[, 1]);

%o {for(i=3,nn, listput(a,1); j=rad(rad(a[i-1])+rad(a[i-2])); while((rad(j)!=rad(rad(a[i-1])+rad(a[i-2])) || #select(n->n==j, a) > 0),j+= rad(rad(a[i-1])+rad(a[i-2])));  a[i]=j; ); print(a);} \\ _Gerry Martens_, Nov 10 2022

%Y Cf. A007947, A005117, A121369, A354184.

%K nonn

%O 1,2

%A _David James Sycamore_, Nov 08 2022

%E More terms from _David A. Corneth_, Nov 09 2022