A358093 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).
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, 46, 61, 107, 42, 149, 191, 170, 19, 21, 20, 961, 41, 18, 47, 53, 40, 63, 29791, 26, 57, 83, 70, 51, 121, 62, 73, 45, 22, 1369, 59, 24, 65, 71, 34, 105, 139, 122, 87, 209, 74
Offset: 1
Keywords
Examples
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. To find a(36), we have 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. [Corrected by _N. J. A. Sloane_, Mar 24 2025] a(43) = 29791 (31^3) because rad(40) + rad(63) = 31; we already have 31 and 31^2.
Links
- Michael De Vlieger, Table of n, a(n) for n = 1..1500
- Michael De Vlieger, Log log scatterplot of a(n), n = 1..1500, highlighting primes in red, other prime powers in gold, and other numbers in blue.
Programs
-
Mathematica
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 *)
-
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 } rad(n) = factorback(factorint(n)[, 1]); \\ David A. Corneth, Nov 09 2022
-
PARI
nn=300; a=List([1,2]); rad(n)= factorback(factorint(n)[, 1]); {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
Extensions
More terms from David A. Corneth, Nov 09 2022
Comments