picoCTF Sequences writeup

picoCTF2022のsequencesが難しかったのでwritteupを書きます。
なお、思い出しながら復習しつつ書いているので実際のコードなどは解いた時と違う可能性があります。

# Sequences
>Description
I wrote this linear recurrence function, can you figure out how to make it run fast enough and get the flag?
Download the code here sequences.py
Note that even an efficient solution might take several seconds to run. If your solution is taking several minutes, then you may need to reconsider your approach.
 添付ファイル
　　sequence.py

***
 sequence.py
```
import math
import hashlib
import sys
from tqdm import tqdm
import functools
ITERS = int(2e7)
VERIF_KEY = "96cc5f3b460732b442814fd33cf8537c"
ENCRYPTED_FLAG = bytes.fromhex("42cbbce1487b443de1acf4834baed794f4bbd0dfe2d6046e248ff7962b")

# This will overflow the stack, it will need to be significantly optimized in order to get the answer :)
@functools.cache
def m_func(i):
    if i == 0: return 1
    if i == 1: return 2
    if i == 2: return 3
    if i == 3: return 4

    return 55692*m_func(i-4) - 9549*m_func(i-3) + 301*m_func(i-2) + 21*m_func(i-1)


# Decrypt the flag
def decrypt_flag(sol):
    sol = sol % (10**10000)
    sol = str(sol)
    sol_md5 = hashlib.md5(sol.encode()).hexdigest()

    if sol_md5 != VERIF_KEY:
        print("Incorrect solution")
        sys.exit(1)

    key = hashlib.sha256(sol.encode()).digest()
    flag = bytearray([char ^ key[i] for i, char in enumerate(ENCRYPTED_FLAG)]).decode()

    print(flag)

if __name__ == "__main__":
    sol = m_func(ITERS)
    decrypt_flag(sol)
```
  decrypt_flag関数によるとENCRYPTED_FLAGはsolのsha256ハッシュとのXORのようなのでsolを求めればよさそうです。solはm_func(ITER)によって生成されます。ITERはイテレーターとか繰り返し回数とか呼ばれるやつで、ここでは2e7回、つまり2×10^7回繰り返されることが変数名から予想できます。実際、m_funcは再帰的に呼び出されています。この回数は計算能力で殴っても無駄だよってことなので、他の方法を考えます。
　m_funcは数列の漸化式になっていて、数学的に言うと

```math
a_n=21a_{n-1}+301a_{n-2}-9549a_{n-3}+55692a_{n-4}
```

　となっているため、この一般項を求めればよさそう。
　一般項を求めるには4次の特性方程式

```math
x^4=21x^3+301x^2-9549x+55692
```

の解が四つの解を持てばよい。
(想像だけど、フィボナッチ数列みたいな無理数でもいいと思う。ただ重解を持つ場合は厳しいと思います。詳しく正確な話が知りたい場合は各種数学のサイトを参考にしてください。)
の解を求める必要がある。
　とりあえず普段使っている[pari/gp](https://pari.math.u-bordeaux.fr/download.html)を使って因数分解してみましょう。


```
gp > f(x)=x^4-21*x^3-301*x^2+9549*x-55692
%1 = (x)->x^4-21*x^3-301*x^2+9549*x-55692
gp > factor(f(x))
%2 =
[x - 17 1]
[x - 13 1]
[x - 12 1]
[x + 21 1]
```
ということで
```math
f(x)=(x-17)*(x-13)*(x-12)*(x+21)\\
```
であることが分かりました。
ここで、x=17を例にして一つ項数を減らしてみます。

```math
f(x)/(x-17)= x^3 - 4x^2 - 369x + 3276\\
\hspace{60pt} a_n=21a_{n-1}+301a_{n-2}-9549a_{n-3}+55692a_{n-4}\\
⇔a_n -4a_{n-1} -369a_{n-2} +3276a_{n-3} \\
\hspace{70pt} =17\big(a_{n-1} -4a_{n-2} -369a_{n-3} +3276a_{n-4} \big)\\
⇔a_n -4a_{n-1} -369a_{n-2} +3276a_{n-3} \\
\hspace{60pt} =17^{n-4}\big(a_4 -4a_3 -369a_2 +3276a_1 \big)\\
=17^{n-4}\big(2350)\\
```
このときの右辺がnに依存するので、同様の手法で項数を減らすことはできません。
数学一般に式を見るときのポイントは、文字数、次数、対称性です。
次数は１次、対称性はないのでこの式が解けない理由は文字数(未知数)が四つもあるためです。
（$nは2*10^7なので定数。a_n,a_{n-1},a_{n-2},a_{n-3}$は未知数）
未知数を求めるためには同じ数の(同値でない)式が必要です。
今回は最初の特性方程式が相違なる４解をもつので、同様のことをあと３回やれば式が四つ得られ$a_n$を求めることができそうです。
そう思って計算していたのですが、手計算にやるのは意外と難しくノートを数十ページ使ってしまったので、簡単な方法を調べます。どうやらsagemathを使えば自動的に説いてくれるそうです。
[Sagemath](https://doc.sagemath.org/html/en/installation/index.html)は格子暗号の問題や、RSAの秘密鍵の一部が欠落している場合の問題(Coppersmithの定理)など、かなり難しい問題を解くときに使うツールというイメージがあり、自分にはまだ敷居が高いと思っていたのですが初挑戦をすることにしました。ちなみに[オンライン版](https://sagecell.sagemath.org/)もありますが、計算に時間がかかる場合は答えが返ってこない場合があり、少なくとも私の環境ではインストール必須でした。
Sagemathは愛用のpari/gpやpythonなども含んでいるので計算結果のSha256を求めるところもスムーズにできそうです。
[ここ](https://qiita.com/taiyaki8926/items/430da4fb7d41c368b9b1)を参考にして連立方程式を解いてみます。ついでに初項の計算もSagemathにしてもらいましょう
```
sage: import hashlib
....: n=int(2e7+1)
....: v1=vector([4,3,2,1])
....: sol = Matrix([[1, -4,-369, 3276],
....:               [1, -8,-405, 4284],
....:               [1, -9,-409, 4641],
....:               [1,-42, 581,-2652]])
....: c =sol*v1
....: v = vector([17^(n-4)*c[0],13^(n-4)*c[1],12^(n-4)*c[2],(-21)^(n-4)*c[3]])
....: keys = sol.solve_right(v)
....: key =keys[0]
....: key =key %(10^10000)
....: key =str(key)
....: m=hashlib.md5()
....: m.update(key.encode())
....: print(m.hexdigest())
96cc5f3b460732b442814fd33cf8537c
```
良さそうです。
求めるべきは$a_{n+1}$なのに注意てITERに1を足しておきます。


```
sage: import hashlib
....: ENCRYPTED_FLAG = bytes.fromhex("42cbbce1487b443de1acf4834baed794f4bbd0dfe2
....: d6046e248ff7962b")
....: n=int(2e7+1)
....: v1=vector([4,3,2,1])
....: sol = Matrix([[1, -4,-369, 3276],
....:               [1, -8,-405, 4284],
....:               [1, -9,-409, 4641],
....:               [1,-42, 581,-2652]])
....: c =sol*v1
....: v = vector([17^(n-4)*c[0],13^(n-4)*c[1],12^(n-4)*c[2],(-21)^(n-4)*c[3]])
....: keys = sol.solve_right(v)
....: key =keys[0]
....: key =key %(10^10000)
....: key =str(key)
....: s=hashlib.sha256()
....: s.update(key.encode())
....: xorkey=s.digest()
....: flag = bytearray([char ^^ xorkey[i] for i, char in enumerate(ENCRYPTED_FLA
....: G)]).decode()
....: print(flag)
picoCTF{b1g_numb3rs_689693c6}
```
解けました。
Sagemath　使いやすいですね。