A354184 - OEIS

%I #27 Aug 12 2022 18:59:16

%S 1,1,2,3,5,8,7,9,10,13,23,36,29,35,64,37,39,76,77,115,192,121,17,28,

%T 31,15,16,17,19,36,25,11,16,13,15,28,29,43,72,49,13,20,23,33,56,47,61,

%U 108,67,73,140,143,213,356,391,569,960,599,629,1228,1243,1857,3100,1867

%N a(1) = a(2) = 1, a(n) = (A007947(31*a(n-1)) + A007947(31*a(n-2)))/31 for n >= 3, i.e., 31*a(n) is the largest squarefree divisor of 31*a(n-1) plus the largest squarefree divisor of 31*a(n-2).

%C After the first 5 terms, the sequence values repeat periodically with a cycle length of 207. The maximum value of a(n) is 1142300, whose first occurrence appears at n = 111.

%H Michael De Vlieger, <a href="/A354184/b354184.txt">Table of n, a(n) for n = 1..2075</a>

%H Augusto Santi, <a href="https://math.stackexchange.com/questions/4452873/">Periodic sequences of integers generated by a(n+1) = rad(a(n)) + rad(a(n-1))</a>, Mathematics Stack Exchange.

%F For n >= 6, a(207 + n) = a(n).

%e 31*2 is the largest squarefree divisor of 31*a(6) = 31*8. 31*7 is the largest squarefree divisor of 31*a(7) = 31*7. So a(8) = (31*2 + 31*7)/31 = 9.

%t Nest[Append[#, (Times @@ FactorInteger[31 #[[-1]]][[All, 1]] + Times @@ FactorInteger[31 #[[-2]]][[All, 1]])/31] &, {1, 1}, 62] (* _Michael De Vlieger_, Jul 18 2022 *)

%o (Python)

%o from sympy import primefactors

%o def rad(num):

%o     primes = primefactors(num)

%o     value = 1

%o     for p in primes:

%o         value *= p

%o     return value

%o a = [1, 1]

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

%o     a += [(rad(31*a[n-1]) + rad(31*a[n-2])) // 31]

%o (PARI) rad(n) = factorback(factorint(n)[, 1]); \\ A007947

%o lista(nn) = my(va = vector(nn)); va[1] = 1; va[2] = 1; for (n=3, nn, va[n] = (rad(31*va[n-2]) + rad(31*va[n-1]))/31;); va; \\ _Michel Marcus_, May 21 2022

%Y Cf. A007947, A121369.

%K nonn,look

%O 1,3

%A _Augusto Santi_, May 18 2022