Rating: 5.0
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:
```python
from 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**-1
print(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 gcd
R=f*R
print(f"f(n) = {L} * {D}**n * {R} // {f}")
n=symbols("n")
print(f"f(n) = ({(L*D**n*R)[0]}) // {f}")
```
Output:
```text
M = 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]]) // 42636
f(n) = (1612*(-21)**n + 981920*12**n - 1082829*13**n + 141933*17**n) // 42636
```
Therefore, in python the function is:
```python
f=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
# OPTIMIZED
from gmpy2 import mpz
f=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.
