login

Year-end appeal: Please make a donation to the OEIS Foundation to support ongoing development and maintenance of the OEIS. We are now in our 61st year, we have over 378,000 sequences, and we’ve reached 11,000 citations (which often say “discovered thanks to the OEIS”).

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