A354184 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).

1, 1, 2, 3, 5, 8, 7, 9, 10, 13, 23, 36, 29, 35, 64, 37, 39, 76, 77, 115, 192, 121, 17, 28, 31, 15, 16, 17, 19, 36, 25, 11, 16, 13, 15, 28, 29, 43, 72, 49, 13, 20, 23, 33, 56, 47, 61, 108, 67, 73, 140, 143, 213, 356, 391, 569, 960, 599, 629, 1228, 1243, 1857, 3100, 1867

Offset: 1

Author

Augusto Santi, May 18 2022

Keywords

• nonn
• look

Comments

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.

Examples

			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.
		

Links

• Michael De Vlieger, Table of n, a(n) for n = 1..2075
• Augusto Santi, Periodic sequences of integers generated by a(n+1) = rad(a(n)) + rad(a(n-1)), Mathematics Stack Exchange.

Crossrefs

Cf. A007947, A121369.

Programs

Mathematica

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

PARI

rad(n) = factorback(factorint(n)[, 1]); \\ A007947
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

Python

from sympy import primefactors
def rad(num):
    primes = primefactors(num)
    value = 1
    for p in primes:
        value *= p
    return value
a = [1, 1]
for n in range(2, 1000):
    a += [(rad(31*a[n-1]) + rad(31*a[n-2])) // 31]

Formula

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