Interesting Crypto LCG problem
from CorCTF 2022
tadpole
没什么复杂的代数变换
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
|
from Crypto.Util.number import bytes_to_long, isPrime
from secrets import randbelow
p = bytes_to_long(open("flag.txt", "rb").read())
assert isPrime(p)
a = randbelow(p)
b = randbelow(p)
def f(s):
return (a * s + b) % p
print("a = ", a)
print("b = ", b)
print("f(31337) = ", f(31337))
print("f(f(31337)) = ", f(f(31337)))
|
exp:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
from Crypto.Util.number import long_to_bytes
a = 7904681699700731398014734140051852539595806699214201704996640156917030632322659247608208994194840235514587046537148300460058962186080655943804500265088604049870276334033409850015651340974377752209566343260236095126079946537115705967909011471361527517536608234561184232228641232031445095605905800675590040729
b = 16276123569406561065481657801212560821090379741833362117064628294630146690975007397274564762071994252430611109538448562330994891595998956302505598671868738461167036849263008183930906881997588494441620076078667417828837239330797541019054284027314592321358909551790371565447129285494856611848340083448507929914
z1 = 52926479498929750044944450970022719277159248911867759992013481774911823190312079157541825423250020665153531167070545276398175787563829542933394906173782217836783565154742242903537987641141610732290449825336292689379131350316072955262065808081711030055841841406454441280215520187695501682433223390854051207100
z2 = 65547980822717919074991147621216627925232640728803041128894527143789172030203362875900831296779973655308791371486165705460914922484808659375299900737148358509883361622225046840011907835671004704947767016613458301891561318029714351016012481309583866288472491239769813776978841785764693181622804797533665463949
s = 31337
f1 = z1 - (a * s + b)
f2 = z2 - (a * (a * s + b) + b)
p = gcd(f1, f2)
print(long_to_bytes(int(p)).decode())
# corctf{1n_m4th3m4t1c5,_th3_3ucl1d14n_4lg0r1thm_1s_4n_3ff1c13nt_m3th0d_f0r_c0mput1ng_th3_GCD_0f_tw0_1nt3g3rs}
|
Luckyguess
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
|
#!/usr/local/bin/python
from random import getrandbits
p = 2**521 - 1
a = getrandbits(521)
b = getrandbits(521)
print("a =", a)
print("b =", b)
try:
x = int(input("enter your starting point: "))
y = int(input("alright, what's your guess? "))
except:
print("?")
exit(-1)
r = getrandbits(20)
for _ in range(r):
x = (x * a + b) % p
if x == y:
print("wow, you are truly psychic! here, have a flag:", open("flag.txt").read())
else:
print("sorry, you are not a true psychic... better luck next time")
|
求解:
形式上来看似乎是一个LCG相关问题
$$
x_{n+1} = x_{n}a + b \pmod{p}
$$
问题在于r = getrandbits(20)
似乎无法预测,然而可以从数列的角度简单理解一下这道题。
这个oracle ,输入初始值$x_{0}$,然后预测$x_{r}$,预测成功并将我们设计得到的$x_{r}$(代码中设计为y
)输入。
关键变量r
也未知,那么可以考虑不动点解法。
先看一组式子:
$$
\begin{cases}
x_{0} &= x_{0} \\
x_{1} &= ax_{0}+b \\
x_{2} &= a^{2}x_{0}+b(a+1) \\
x_{3} &= a^{3}x_{0}+b(a^2+a+1)
\end{cases}
$$
$$
x_{n} = a^{n}x_{0}+ b\left( \sum^{n-1}_{i=0}a^{i} \right)
$$
注意到
$$
\frac{a^{n}-1}{a-1}= \left( \sum^{n-1}_{i=0} a^{i} \right)
$$
,
实际上可以写出一个较为简洁的关于第$n$项和首项的公式:
$$
x_{n} = a^{n}x_{0}+ b\frac{a^{n}-1}{a-1}
$$
这就符合我们控制的两个输入和输出。
首先不考虑模运算的存在,一个较好的情况是$x_{n}=x_{0},n=r$,$r$未知,那么考虑消除$r$对我们的影响:
$$\begin{aligned}
\frac{x_{r}}{x_{0}}&=a^{r}+\frac{b}{x_{0}} \frac{a^{r}-1}{a-1}\\
1&=a^{r}+\frac{b}{x_{0}} \frac{a^{r}-1}{a-1}
\end{aligned}$$
OK,到了这里我们发现等式右边由$a^{n}$和一个系数$\frac{b}{x_{0}}$控制的$a^{n}$构成,那么一个平凡的想法就是考虑两个系数互为相反数,那么尝试一下我们刚好发现等式恒成立。
令$\frac{b}{x_{0}} \frac{1}{a-1}=-1$
$$x_{0}= - \frac{b}{a-1}$$
代入恰好得:
$$1 = a^{r}-(a^{r}-1)=1$$
那么此时,只要已知$a,b$令$x_{0}=- \frac{b}{a-1}$,这个递推数列就变成了常数数列,再考虑模运算自然也是无伤大雅。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
from pwn import *
conn = remote('be.ax', 31800)
# conn = process("name.py")
p = 2**521 - 1
a = int(conn.recvline().decode().strip().split('a = ')[1])
b = int(conn.recvline().decode().strip().split('b = ')[1])
x = y = -b * pow(a - 1, -1, p) % p
conn.sendlineafter(b'starting point: ', str(x).encode())
conn.sendlineafter(b'guess? ', str(y).encode())
print(conn.recvline().decode())
# corctf{r34l_psych1c5_d0nt_n33d_f1x3d_p01nt5_t0_tr1ck_th15_lcg!}
|
妙!
但是推导到这里,似乎有些凑巧。
先看下一题:
Exchange
Task
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
|
from Crypto.Util.number import *
from Crypto.Cipher import AES
from Crypto.Util.Padding import pad
from hashlib import sha256
from secrets import randbelow
p = 142031099029600410074857132245225995042133907174773113428619183542435280521982827908693709967174895346639746117298434598064909317599742674575275028013832939859778024440938714958561951083471842387497181706195805000375824824688304388119038321175358608957437054475286727321806430701729130544065757189542110211847
a = randbelow(p)
b = randbelow(p)
s = randbelow(p)
print("p =", p)
print("a =", a)
print("b =", b)
print("s =", s)
a_priv = randbelow(p)
b_priv = randbelow(p)
def f(s):
return (a * s + b) % p
def mult(s, n):
for _ in range(n):
s = f(s)
return s
A = mult(s, a_priv)
B = mult(s, b_priv)
print("A =", A)
print("B =", B)
shared = mult(A, b_priv)
assert mult(B, a_priv) == shared
flag = open("flag.txt", "rb").read()
key = sha256(long_to_bytes(shared)).digest()[:16]
iv = long_to_bytes(randint(0, 2**128))
cipher = AES.new(key, AES.MODE_CBC, iv=iv)
print(iv.hex() + cipher.encrypt(pad(flag, 16)).hex())
|
利用第二题推导的通项公式,
$$\begin{aligned}
A &= f_{n_{a}}(s)= a^{n_{a}}s+ b \frac{a^{n_{a}}-1}{a-1} \pmod{p}\\
A(a-1)+b&=a^{n_{a}}(as-s+b) \pmod{p}\\
a^{n_{a}} &= \frac{A(a-1)+b}{as-s+b} \pmod{p}
\end{aligned}$$
从上式的推导可以看出,$a^{n_{a}},a^{n_{b}}$可以由已知量导出。
这就回到了DH的一般思考场景。
验证得,$p-1$光滑,直接求解DLP,计算
$$
\begin{aligned}
n_{a} &= DiscreteLog(a^{n_{a}},a,p)\\
shared &= f_{n_{a}}(B)\\
&= a^{n_{a}}B+b \frac{a^{n_{a}}-1}{a-1} \pmod{p}
\end{aligned}
$$
唯一未知量已可求解,获取了$shared$
可以求解。
exp:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
|
sage: from Crypto.Util.number import long_to_bytes
....: from Crypto.Cipher import AES
....: from Crypto.Util.Padding import unpad
....: from hashlib import sha256
....:
....: p = 1420310990296004100748571322452259950421339071747731134286191835424352805219828279086937099671748953466397461
....: 17298434598064909317599742674575275028013832939859778024440938714958561951083471842387497181706195805000375824824
....: 688304388119038321175358608957437054475286727321806430701729130544065757189542110211847
....: a = 1180906598237265321184570154603935013535512571819012348308688052993667257580121658456389778783222827629290215
....: 70278435511082796994178870962500440332899721398426189888618654464380851733007647761349698218193871563040337609238
....: 025971961729401986114391957513108804134147523112841191971447906617102015540889276702905
....: b = 5795014987100615243467302014637519655589220562695967625172441001618493582571250812112330936022277755982709396
....: 54689652681477200276478424926550717060636693281351272022500409354148364163603509242184627980038782665632058932676
....: 35176851677889275076622582116735064397099811275094311855310291134721254402338711815917
....: s = 3570158135111160465491334886700707833940269177041036813362503042720279105776685310351097408959241134406576995
....: 73708026173784951618374426701578277686774118710424015000713663174396814612714838808580074695024533617060019734419
....: 02698612564888892738986839322028935932565866492285930239231621460094395437739108335763
....: A = 2705569950255528261367920540242672730435988633782267523285646370856059877266600466366005252832869228207716559
....: 02594950903882166292400533970414295870526111331638869384711648295375897115982531152701610900861800015012271649251
....: 99272064309777701514693535680247097233110602308486009083412543129797852747444605837628
....: B = 1321783200371127370097264683674718982421959235681582348717736070054240011526943389939787036890301472158431250
....: 95282272730052868843423659165019475476788785426513627877574198334376818205173785102362137159225281640301442638067
....: 549414775820844039938433118586793458501467811405967773962568614238426424346683176754273
....: ct = bytes.fromhex('e0364f9f55fc27fc46f3ab1dc9db48fa482eae28750eaba12f4f76091b099b01fdb64212f66caa6f366934c3b9929
....: bad37997b3f9d071ce3c74d3e36acb26d6efc9caa2508ed023828583a236400d64e')
....:
sage: factor(p-1)
2 * 593603 * 635603 * 904901 * 910981 * 994141 * 1270897 * 1432181 * 1612759 * 1644901 * 1679267 * 1766573 * 1777031 * 1785209 * 2014487 * 2088217 * 2733233 * 3101729 * 5796797 * 7243927 * 7731943 * 9403621 * 15787589 * 23697101 * 24045389 * 29428913 * 31794713 * 34572653 * 37082147 * 62263273 * 76785581 * 85536359 * 112094921 * 125101001 * 132504719 * 133775189 * 155631871 * 172174381 * 180502669 * 251177879 * 347924039 * 428295697 * 466014559 * 501589381
sage: # p-1 smooth, dlp go
sage: F = GF(p)
sage: ax = F(b+ (a-1)*A) / F(a*s - s + b)
sage: x = discrete_log(F(ax),F(a))
sage: x
63497966771228335993935218724355716676359926182967975463093060105894867914802051403524263838451533117998140812842659902100945139428076742018151358770711075499718955791059569760306592803008376586431708302909649602374612770031446609941268056564386982932718212036534269872682126736591975854896102855270700008372
sage: shared = F(a)^x * B + b * F(F(a)^x - 1) / F(a - 1)
....:
....: key = sha256(long_to_bytes(int(shared))).digest()[:16]
....: iv = ct[:16]
....: cipher = AES.new(key, AES.MODE_CBC, iv=iv)
....: flag = unpad(cipher.decrypt(ct[16:]), 16)
....: print(flag.decode())
....:
corctf{th1s_lcg_3xch4ng3_1s_4_l1ttl3_1ns3cur3_f0r_n0w}
|
leapfrog
本题是对第一题的扩展,本题的任务是根据一些函数输出恢复$a,b,p$.
针对不连续的输出值,我们有以下处理方式
#task
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
|
from Crypto.Util.number import long_to_bytes, getPrime
from Crypto.Cipher import AES
from Crypto.Util.Padding import pad
from hashlib import sha256
from secrets import randbelow
from random import sample
p = getPrime(256)
a = randbelow(p)
b = randbelow(p)
s = randbelow(p)
def f(s):
return (a * s + b) % p
jumps = sample(range(3, 25), 12)
output = [s]
for jump in jumps:
for _ in range(jump):
s = f(s)
output.append(s)
print(jumps)
print(output)
flag = open("flag.txt", "rb").read()
key = sha256(b"".join([long_to_bytes(x) for x in [a, b, p]])).digest()[:16]
iv = long_to_bytes(randbelow(2**128))
cipher = AES.new(key, AES.MODE_CBC, iv=iv)
print(iv.hex() + cipher.encrypt(pad(flag, 16)).hex())
|
1
2
3
4
|
jumps = [5, 3, 23, 13, 24, 6, 10, 9, 7, 4, 19, 16]
output = [26242498579536691811055981149948736081413123636643477706015419836101346754443, 30320412755241177141099565765265147075632060183801443609889236855980299685595, 65684356693401962957802832810273549345608027337432965824937963429120291339333, 15025547765549333168957368149177848577882555487889680742466312084547650972663, 46764069432060214735440855620792051531943268335710103593983788232446614161424, 71575544531523096893697176151110271985899529970263634996534766185719951232899, 8149547548198503668415702507621754973088994278880874813606458793607866713778, 12081871161483608517505346339140143493132928051760353815508503241747142024697, 65627056932006241674763356339068429188278123434638526706264676467885955099667, 23413741607307309476964696379608864503970503243566103692132654387385869400762, 56014408298982744092873649879675961526790332954773022900206888891912862484806, 77000766146189604405769394813422399327596415228762086351262010618717119973525, 14589246063765426640159853561271509992635998018136452450026806673980229327448]
iv_c = 05ac5b17c67bcfbf5c43fa9d319cfc4c62ee1ce1ab2130846f776e783e5797ac1c02a34045e4130f3b8111e57397df344bd0e14f3df4f1a822c43c7a89fd4113f9a7702b0b0e0b0473a2cbac25e1dd9c
|
为了便于表示,我们引入$J_{i}$表示jumps
中的第$i$个数,$S_{i}$表示jumps
中的前$n$项和,显然地,我们得到:
$J_{1}=5,J_{2}=3,\dots;S_{0}=0,S_{1}=5,S_{2}=3\dots$
根据
1
2
3
4
5
6
7
8
9
|
def f(s):
return (a * s + b) % p
jumps = sample(range(3, 25), 12)
output = [s]
for jump in jumps:
for _ in range(jump):
s = f(s)
output.append(s)
|
根据前期的推导:
$$f_{n}(x_{0})=a^{n}x_{0}+b \frac{a^{n}-1}{a-1} \pmod{p}$$
可以得到output
中的表达结果为(如果记output列表为${ O_{i} }$:
$$
O_{i}=f_{S_{i}}(s)= a^{S_{i}}s+ b \frac{a^{S_{i}}-1}{a-1} \pmod{p}
$$
那么对于12=len(jumps)
个输出值和指数,我们得到了一簇$l=12$个方程。
在$\mod{p}$的情况下,常见的思路就是加减消元。
不难观察到,
$$a^{69} = \frac{O_{9}-O_{6}}{O_{3}-O_{1}} \pmod{p} \tag{1}$$
如果能够找到另外一个$a^{k}$,就给了我们进一步作差的思路:
如果能够找到两个式子$f_{1},f_{2}$满足:
$$\begin{cases}
f_{1} = 0 \\
f_{2} = 0
\end{cases} \pmod{p}$$
那么可以获得$k_{1}p=(f_{1},f_{2})$,那么此时的$k_{1}$相对小,可以通过可以接受的枚举,恢复出最终模数$p$.
对于$(1)$式,这样的一组$a^k,i_{1},i_{2},i_{3},i_{4}$
对于$l=12$的情况,考虑枚举,来获得我们想要的
$$a^{k}= \frac{O_{i_{1}}-O_{i_{2}}}{O_{i_{3}}-O_{i_{3}}} \pmod{p}$$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
|
sage: jumps = [5, 3, 23, 13, 24, 6, 10, 9, 7, 4, 19, 16]
....:
....: def f(i):
....: return a^i * s + b * sum(a^j for j in range(i))
....:
....: P.<a, s, b> = PolynomialRing(QQ)
....: S = [0] + [sum(jumps[:i]) for i in range(1, len(jumps) + 1)]
....: Z = [f(s) for s in S]
....:
....: good = []
....: for i1, i2, i3, i4 in Combinations(range(len(Z)), 4):
....: f_ = (Z[i4] - Z[i3])/(Z[i2] - Z[i1])
....: if f_.denominator() == 1 and len(P(f_).coefficients()) == 1:
....: good.append((f_, i1, i2, i3, i4))
....:
....: print(good)
[(a^69, 1, 3, 6, 9), (a^79, 1, 4, 7, 11), (a^95, 1, 4, 9, 12), (a^92, 2, 3, 9, 11), (a^60, 2, 4, 5, 10), (a^49, 4, 6, 8, 11), (a^55, 5, 7, 11, 12), (a^30, 6, 8, 10, 11), (a^39, 7, 9, 11, 12)]
|
欣喜地,观察到
$$\begin{cases}
a^{69}=a^{39}a^{30} \\
a^{60}=a^{30}a^{30} \\
a^{79}=a^{30}a^{49}
\end{cases}$$
那么符合条件的$f_{1},f_{2},f_{3}$已找到
恢复$p$:
1
2
3
4
5
6
7
|
#sage
sage: f1 = (O[9] - O[6]) * (O[8] - O[6]) * (O[9] - O[7]) - (O[11] - O[10]) * (O[12] - O[11]) * (O[3] - O[1])
....: f2 = (O[10] - O[5]) * (O[8] - O[6])^2 - (O[11] - O[10]) ^ 2 * (O[4] - O[2])
sage: p = ZZ(gcd(f1, f2))
sage: factor(p)
2^2 * 3^3 * 82854189798454303629690440183768928075006051188051668052925581650741089039941
sage: p = 82854189798454303629690440183768928075006051188051668052925581650741089039941
|
那么对于恢复了$p$之后的求解$a,b,s$过程无非是对于一簇$12$个方程在有限域$GF(p)$上构成的方程组求解其中三个未知量的过程,这不是一件困难的事情。
这里可以应用结式(resultant)求解,颇为优雅。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
|
#sage
from Crypto.Util.number import long_to_bytes
from Crypto.Cipher import AES
from Crypto.Util.Padding import unpad
from hashlib import sha256
from sage.matrix.matrix2 import Matrix
def resultant(f1, f2, var):
return Matrix.determinant(f1.sylvester_matrix(f2, var))
jumps = [5, 3, 23, 13, 24, 6, 10, 9, 7, 4, 19, 16]
outputs =
ct = bytes.fromhex('05ac5b17c67bcfbf5c43fa9d319cfc4c62ee1ce1ab2130846f776e783e5797ac1c02a34045e4130f3b8111e57397df344bd0e14f3df4f1a822c43c7a89fd4113f9a7702b0b0e0b0473a2cbac25e1dd9c')
def f(i):
return a^i * s + b * sum(a^j for j in range(i))
P.<a, s, b> = PolynomialRing(QQ)
S = [0] + [sum(jumps[:i]) for i in range(1, len(jumps) + 1)]
O = outputs
f1 = (O[9] - O[6]) * (O[8] - O[6]) * (O[9] - O[7]) - (O[11] - O[10]) * (O[12] - O[11]) * (O[3] - O[1])
f2 = (O[10] - O[5]) * (O[8] - O[6])^2 - (O[11] - O[10]) ^ 2 * (O[4] - O[2])
p = ZZ(gcd(f1, f2))
#for d in range(2, 0x1000):#The effect of k1 can be eliminated by enumeration, and p' can also be factorized directly under suitable conditions
# while p % d == 0:
# p //= d
p = max(p)[0]
assert p.is_prime()
s = outputs[0]
P.<a, b> = PolynomialRing(GF(p))
f1 = f(S[1]) - outputs[1]
f2 = f(S[2]) - outputs[2]
h1 = resultant(f1, f2, b)
h2 = resultant(f1, f2, a)
a = h1.univariate_polynomial().roots()[0][0]
b = h2.univariate_polynomial().roots()[0][0]
key = sha256(b"".join([long_to_bytes(int(x)) for x in [a, b, p]])).digest()[:16]
cipher = AES.new(key, AES.MODE_CBC, iv=ct[:16])
flag = unpad(cipher.decrypt(ct[16:]), 16)
print(flag.decode())
# corctf{:msfrog:_is_pr0ud_0f_y0ur_l34pfr0gg1ng_4b1lit135}
|
LCG
Linear Congruential Generator
$$X_{n+1} = aX_{n}+b \pmod{p}$$
$A,B,P$是产生器设定的常数。
LCG的周期最大为$p$,但大部分情况下均小于$p$,要令$LCG$达到最大周期,符合以下条件:
- $(b,p)=1$
- $M$的所有质因数都能整除$a-1$
- 若$p$是$4$的倍数,$a-1$也是
- $a,b,X_{0}$均小于$M$
- $a,b$是正整数
那么针对上述一系列题目的一些纵深思考:
递推数列与递推函数
几乎类似的我们在伪随机和数列中都遇见过迭代函数:
$$a_{n+1}=f(a_{n})$$
不难发现,求解递推数列和迭代函数本质并无差别。
记$n$次迭代函数
$$f_{n}(x)= \underbrace{f(f(f( \dots f}_{n}(x))))$$
并补充定义$f_{0}(x)=x$
那么对于递推数列$a_{n+1}=f(a_{n})$也可以通过求解迭代函数的方法求解数列$a_{n}=f_{n}(a_{0})$
函数的不动点
假设有非空集合$A$以及集合映射$f:A \rightarrow A$,我们称$x_{0} \in A$是函数$f$的不动点,如果$x_{0}$满足:
$$f(x_{0})=x_{0}$$
- 与递推式形式一致的函数也被叫做递推关系式的特征函数(characteristic function)或背景函数。
- 递推式特征函数的不动点也叫做递推式的不动点。不动点在对应函数图象上的对应点也叫做函数图象的不动点
$f$的不动点的集合称为不动点集。
特殊的,如果函数$g$为可逆函数:
- $(x_{0},y_{0})$为函数$g$的不动点,那么$(y_{0},x_{0})$为函数$g^{-1}$的不动点
- 函数不动点坐标具有$(x_{0},x_{0})$的形式,是函数图像与$y=x$的交点.
- 将$g$与$g^{-1}$的图像绘制在同一坐标系中,两图像交点均为不动点
不动点的概念源于拓扑。
但是他的用途广泛,一直是一个有趣的数学概念。
比如巴拿赫不动点定理,又称压缩映射定理,这个定理指出,如果一张地图的内容包含地图本身所处的位置,那么地图上一定有一个点和它的实际位置重合,也就是“从现实位置到地图位置”这一映射的不动点。
比如数值分析中的不动点迭代法,可以快速求解一个函数的不动点的近似值,进而也就能求解方程的根的近似值。
不动点方法
一阶线性递推数列
对于形如$a_{n+1}=pa_{n}+q$,其中$p,q$为常数,$p \ne 1$称为一阶常系递推数列。
给出一个简化过程;
先设$a_{n+1}=ka_n+d$
再假设我们期待的构造好的等比数列形式$a_{n+1}-p=k(a_n-p)$
用后一式的两端同时对前一式的两端作差,可得:
$(a_{n+1}-p)-a_{n+1}=k(a_n-p)-(ka_n+d)$
$$\Rightarrow-p=-kp-d\Rightarrow\quad p=kp+d$$
最后一个式子显然说明所需要的常数$p$就是函数$f(x)=kx+d$的不动点。
从几何上讲,从截距不为零的一次函数(对应于上述非齐次一阶线性递推关系式)变形为正比例函数(对应于构造出来的等比数列满足的关系式)的变换就是一种图象平移操作,而此类函数的不动点提供了同时沿y轴方向和x轴方向进行图象平移的距离参考值。
前面的推到已经证明了这样的数列的通项公式为
$$a_{n}=a_{1}p^{n-1}+ q \frac{1-p^{n-1}}{1-p}$$
一阶常系数分式线性递推关系式的概念及其通项的待定系数解法
一阶常系数分式线性递推关系式,形如:$$a_{n+1}= \frac{Aa_{n}+B}{Ca_{n}+D}(n \in \mathbb{N^{+}},A,B,C,D \in \mathbb{R},AD-BD \ne 0)$$
不动点方法:
设函数$f(x)= \frac{Ax+B}{Cx+D},x,A,B,C,D \in \mathbb{C},AD-BC \ne 0$,其不动点为
$$x=\lambda=f(\lambda)\tag{2}$$
令$x_{0}$为$(2)$的任意一个根,从而$x_{0}$是一个不动点;
$$\begin{aligned}
a_{n+1}-x_{0} &= \frac{Aa_{n}+B}{(Ca_{n}+D)}-x_{0} \\
&=\frac{(A-Cx_{0})(a_{n}-x_{0})}{Ca_{n}+D} \\
\frac{1}{a_{n+1}-x_{0}}&= \frac{Ca_{n}+D}{(A-Cx_{0})(a_{n}-x_{0})} \\
&=\frac{Cx_{0}+D}{A-Cx_{0}} \frac{1}{a_{n}-x_{0}}+ \frac{C}{A-Cx_{0}}
\end{aligned}$$
由上式可知,
$$\set{ \frac{1}{a_{n}-x_{0}} }$$
整体来说,是一个一阶常系数线性递推数列,且包含有取值为常值的非齐次项。
对于关于$\lambda$方程$\lambda=f(\lambda)$,对其解进行分类讨论:
$i.$
方程$(2)$的两个根为重根:$x_{1}=x_{2}= \frac{A-D}{2C}$
特别地,对于$\frac{Cx_{1}+D}{A-Cx_{1}}=1$时,可以继续化简:
$$\frac{1}{a_{n+1}-x_0}=\frac{Cx_0+D}{A-Cx_0}\frac{1}{a_n-x_0}+\frac C{A-Cx_0}=\frac{1}{a_n-x_0}+\frac C{A-Cx_0}$$
即:$\frac{1}{a_{n}-x_{0}}$变为以$\frac{1}{a_{n}-x_{0}}$为首项,$\frac{C}{A-Cx_{0}}$为公差的等差数列,容易知道通项公式:
$$\frac{1}{a_{n}-x_{0}}=\frac{1}{a_{n}-x_{0}}+\frac{2C(n-1)}{A+D}$$
从而容易的推出$a_{n}$的解析式。
$ii.$
如果不动点方程$(2)$有两个相异的根$x_{1},x_{2}$,那么类似地,考虑两侧同时减去不动点的思路:
$$
\begin{cases}
a_{n+1}-x_{1}= \frac{(A-Cx_{1})(a_{n}-x_{1})}{Ca_{n}+D}\\
a_{n+1}-x_{2}= \frac{(A-Cx_{2})(a_{n}-x_{2})}{Ca_{n}+D}
\end{cases}
$$
考虑分母相同,两式作商:
$$\frac{a_{n+1}-x_{1}}{a_{n+1}-x_{2}}=k \cdot \frac{a_{n}-x_{1}}{a_{n}-x_{2}},k=\frac{A-Cx_{1}}{A-Cx_{2}}$$
易知,$\set{\frac{a_{n}-x_{1}}{a_{n}-x_{2}}}$是以$\frac{a_{1}-x_{1}}{a_{1}-x_{2}}$为首项,$k$为公比的等比数列
代换得出$a_{n}$的解析式并不困难。
那么他们从何而来?
分式线性变换的概念与分式线性递推式的不动点解法由来
定义:具有以下形式的函数$f:\mathbb{C} \rightarrow \mathbb{C}$,称为分式线性变换(linear fractional transformation):
$$f(x) = \frac{Ax+B}{Cx+D}\quad(A,B,C,D \in \mathbb{C},AD-BC \ne 0)$$
因为德国数学家奥古斯特·莫比乌斯曾对其进行过大量研究,所以分式线性变换也叫做莫比乌斯变换(Möbius transformation)
因式分解后可把$f(x)$写成:
$$f(x)=\frac{BC-AD}C\frac1{Cx+D}+\frac AC$$
可以验证在实数域上,这是一个可逆变换,具有以下性质:
- $y=f(x)$的逆变换为$x= \frac{-Dy+B}{Cy-A}$。
- 一次函数是一个特殊情况。
- 分式线性变换之间的复合结果仍然是分式线性变换。
- 根据不动点的性质,对于变换及其逆变换,不动点一一对应。
为什么两边同时减去不动点?
按照待定系数法的思路,假定对分式线性变换解析式的$f(x)$两侧同时减去某个合适的常数P再同时取倒数后,能得到形式一致的分母。先
$$
\begin{aligned}&\text{做第一步减法,得:}\\&a_{n+1}-P=\frac{Aa_n+B}{Ca_n+D}-P=\frac{Aa_n+B-PCa_n-DP}{Ca_n+D}=\frac{(A-PC)a_n+(B-DP)}{Ca_n+D}\end{aligned}
$$
要使上面的式子取倒数后的分母满足要求,由于分母取倒数前是位于分子的位置,基于前面的解题实例可知,我们只需要保证上式最左边的量与最右边的分子形式相似即可。
即只需要使$a_{n+1}-P$与$(A-PC)a_n+(B-DP)$的对应系数成相同比例即可。即:
$$\begin{aligned}\frac{1}{A-PC}=\frac{-P}{B-DP}\end{aligned}$$
$$P=\frac{B-DP}{CP-A}$$
显然,这个所需的P值是$f(x)$的反函数的不动点。根据前面总结的性质,可知P也是$f(x)$的不动点。
即对$f(x)$两侧同时减去其不动点后,再取同时倒数,就一定能构造出一次的常系数非齐次递推式。
桥函数
上述推导,无论赛题或概念都在中间利用了大量中间函数。
考虑到分式线性变换性质3.,我们继续探究:
首先由于一次函数属于特殊的分式线性变换,而且分式线性变换是可逆变化,所以自然想到能否将分式线性递推式还原为某种更简单的一次线性递推式的形式。求出与之相关的一次线性递推式的迭代结果后,再利用逆变换反推出原递推式的多次迭代结果。
另一方面,我们已经知道使用不动点法可以简化一次线性递推式的求解,我们希望将不动点也整合到分式线性变换还原为线性递推式的系数配凑过程中。
定义:
使得$f(x)$可以与另一个函数$g(x)$建立下列联系的可逆函数$T(x)$叫做桥函数:
$$f(x)= T^{-1}(g(T(x)))$$
此时,$f(x)$和$g(x)$叫做关于桥函数$T(x)$的一对相似函数,或者说$f(x)$和$g(x)$关于桥函数$T(x)$互为相似函数。
用拓扑学术语来说,桥函数是一种同胚变换(homeomorphism)
它建立了$f(x)$和$g(x)$之间的拓扑共轭(topological conjugacy)
笔者的拓扑记忆也是较为模糊,如有疏漏欢迎勘误。
桥函数是计算函数迭代或构造数列的辅助函数。
通过找到合适的桥函数,将不易得到迭代公式的$f(x)$转化为容易得到迭代公式的形式更简单的$g(x)$, 计算完$g(x)$的$n$次迭代后,再使用桥函数的逆变换(性质保证逆变换一定存在)。应用这种思路把函数迭代问题转化为寻找一个良好的桥函数问题。
那么很容易的,可以得到使用不动点法的动机:
借住不动点,可以探测目标桥函数的一些性质,以便确定桥函数的待定系数。只要确定了待定系数的细节,就可以通过相似变换的方法得到原函数的$n$次迭代式。
举例
假设$f(x)=T^{-1}(g(T(x)))$,且$x_{0}$是$f(x)$的不动点,那么$T(x_{0})$一定也会是$g(x)$的不动点。这说明,如果求出$f(x)$的不动点,并且找到结构合适的$T(x)$,那么$g(x)$的不动点是不言自明的。
通常为了便于求解$g(x)$的$n$次迭代式,可以尝试将$g(x)$取为如下形式:
$$\begin{aligned}
&\bullet g(x) =x+a \\
&\bullet g(x) =ax \\
&\bullet g(x) =ax^2 \\
&\bullet g(x) =ax^3
\end{aligned}$$
对于这几种情况,$g(x)$的不动点为0或$\infty$。此时如果$f(x)$只有唯一的不动点a,可以考虑取$g(x)=x-a$或$g(x)=\frac1{x-a}$;如果$f(x)$有2个相异的不动点a和b,则可以考虑取$g(x)=\frac{x-a}{x-b}$。
这样的桥函数思想是优美的。
另举一例 (from TSG-CTF 2023)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
|
import os
import random
import string
flag = os.getenv("FLAG", "FAKECTF{THIS_IS_FAKE}")
key = [random.randrange(10 ** 4) for _ in flag]
cs = string.printable[:-6]
def r(k):
for _ in range(k):
random.seed(x := random.randrange(20231104, 20231104 * 10))
return x
random.seed(int(input("seed: ")))
print('flag:', ''.join([cs[(cs.index(f) + r(k)) % len(cs)] for f, k in zip(flag, key)]))
|