It is possible to restate the function in the following form (for full details see [https://www.nullhardware.com/reference/hacking-101/picoctf-2022-greatest-hits/sequences/](https://www.nullhardware.com/reference/hacking-101/picoctf-2022-greatest-hits/sequences/)):```                    n                         f(n) = (1 0 0 0) . M  . x                                              o                     where                                              /0     1     0   0 \                          |                  |                          |0     0     1   0 |                      M = |                  |                          |0     0     0   1 |                          |                  |                          \55692 -9549 301 21/                       and                                                 /1\                                           | |                                           |2|                                      x  = | |                                       0   |3|                                           | |                                           \4/                                       Therefore                                                                             n                         /0     1     0   0 \    /1\                   |                  |    | |                   |0     0     1   0 |    |2|f(n) = (1 0 0 0) . |                  |  . | |                   |0     0     0   1 |    |3|                   |                  |    | |                   \55692 -9549 301 21/    \4/
```
You can then use eigendecomposition to factor M and directly derive an equation in terms of n:
```pythonfrom sympy import *M=Matrix([[0,1,0,0],[0,0,1,0],[0,0,0,1],[55692,-9549,301,21]])P,D = M.diagonalize()Pi=P**-1print(f"M = {M}\nP*D*P^-1 = {P*D*Pi}\n")L=Matrix([[1,0,0,0]])*P # pre-multiply by [1,0,0,0]R=Pi*Matrix([[1],[2],[3],[4]]) #post-multiply by [1;2;3;4]f=1/gcd(tuple(R)) # pull out the gcdR=f*Rprint(f"f(n) = {L} * {D}**n * {R} // {f}")n=symbols("n")print(f"f(n) = ({(L*D**n*R)[0]}) // {f}")```
Output:```textM = Matrix([[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [55692, -9549, 301, 21]])P*D*P^-1 = Matrix([[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [55692, -9549, 301, 21]])
f(n) = Matrix([[-1, 1, 1, 1]]) * Matrix([[-21, 0, 0, 0], [0, 12, 0, 0], [0, 0, 13, 0], [0, 0, 0, 17]])**n * Matrix([[-1612], [981920], [-1082829], [141933]]) // 42636f(n) = (1612*(-21)**n + 981920*12**n - 1082829*13**n + 141933*17**n) // 42636```
Therefore, in python the function is:```pythonf=lambda n: (1612*((-21)**int(n)) + 981920*((12)**int(n)) - 1082829*((13)**int(n)) + 141933*((17)**int(n)))//42636```
Which works, but is still slow untill you replace the large integers with a more optimal type (mpz from gmpy2)```python# OPTIMIZEDfrom gmpy2 import mpzf=lambda n: (1612*(mpz(-21)**int(n)) + 981920*(mpz(12)**int(n)) - 1082829*(mpz(13)**int(n)) + 141933*(mpz(17)**int(n)))//42636```
Replacing `m_func` for `f` should give you the answer in approx 1 second. 