login
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).
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,3
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.
LINKS
FORMULA
For n >= 6, a(207 + n) = a(n).
EXAMPLE
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.
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 *)
PROG
(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]
(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
CROSSREFS
Sequence in context: A116918 A271621 A116917 * A121369 A284172 A230445
KEYWORD
nonn,look
AUTHOR
Augusto Santi, May 18 2022
STATUS
approved