login
a(n+1) = 1 + number of previous terms that share a factor > 1 with a(n); a(1) = 2.
1

%I #35 Aug 19 2023 00:13:36

%S 2,2,3,2,4,5,2,6,8,8,9,4,10,12,14,13,2,14,15,8,16,17,2,18,22,20,23,2,

%T 22,23,3,8,24,29,2,26,28,28,29,3,10,32,31,2,32,33,13,4,34,36,42,43,2,

%U 38,39,17,4,40,43,3,15,21,21,22,43,4,43,5,9,19,3,20,48,58,48,60,63,28,52,53,2,51,28

%N a(n+1) = 1 + number of previous terms that share a factor > 1 with a(n); a(1) = 2.

%C There are prominent lines that have more terms, their coefficients are approximately: 0.519, 0.329, 0.689, 0.188, 0.615, ... (see the frequency link). They seem to be distorted prime harmonic lines: 1/2, 1/3, 2/3, 1/5, 3/5, ... from A016035.

%C It appears limsup a(n)/n is approximately 0.83.

%H Rok Cestnik, <a href="/A364934/b364934.txt">Table of n, a(n) for n = 1..1000</a>

%H Rok Cestnik, <a href="/A364934/a364934.pdf">a(n)/n frequency for 3*10^5 terms</a>.

%H Michael De Vlieger, <a href="/A364934/a364934.png">Log log scatterplot of a(n)</a>, n = 1..2^16.

%e [2,*] 1 term shares a factor with 2, so a(2) = 1+1 = 2.

%e [2,2,*] 2 terms share a factor with 2, so a(3) = 1+2 = 3.

%e [2,2,3,*] 1 term shares a factor with 3, so a(4) = 1+1 = 2.

%e [2,2,3,2,*] 3 terms share a factor with 2, so a(5) = 1+3 = 4.

%e [2,2,3,2,4,*] 4 terms share a factor with 4, so a(6) = 1+4 = 5.

%t nn = 120; c[_] := 0; s = {}; a[1] = j = 2; Do[c[j]++; If[c[j] == 1, AppendTo[s, j]]; k = 1 + Total@ Map[c[#] Boole[GCD[#, j] > 1] &, s]; Set[{a[n], j}, {k, k}], {n, 2, nn}]; Array[a, nn] (* _Michael De Vlieger_, Aug 16 2023 *)

%o (Python)

%o from sympy.ntheory import primefactors

%o a=[2]; p = [{2}];

%o for n in range(1,1000):

%o a.append(sum(1 if p[-1].intersection(p[i]) else 0 for i in range(n))+1)

%o p.append(set(primefactors(a[-1])))

%o (Python)

%o from math import gcd, prod

%o from sympy import factorint

%o from itertools import islice

%o from collections import Counter

%o def agen(): # generator of terms

%o a=[2]; f=2; c=Counter([f])

%o while True:

%o yield a[-1]

%o a.append(1 + sum(c[fi] for fi in c if gcd(f,fi)>1))

%o f=prod(factorint(a[-1]).keys())

%o c[f] += 1

%o a=list(islice(agen(), 1000)) # _Michael S. Branicky_, Aug 16 2023

%Y Cf. A016035, A124056.

%K nonn

%O 1,1

%A _Rok Cestnik_, Aug 15 2023