解題區/解題區 - Zerojudge

b443. 我愛 Fibonacci

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution
    1. 4.1. 參考資料
    2. 4.2. 故事

Problem

$$F_n=\begin{cases} n, & n=0,1 \\ F_{n-1}+F_{n-2}, & n \geq 2 \end{cases}$$

求出$F_{2^n} \mod m$ 的結果。

Sample Input

1
2
3
4
3
1 1000000007
2 1000000007
3 1000000007

Sample Output

1
2
3
1
3
21

Solution

一般的矩陣計算，利用 $M^n$ 求出$F_n$，其中

$$M = \begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix}$$

如果是求$F_n$ 時間複雜度 $O(\log n)$，而這一題求的是$F_{2^n}$，時間複雜度 $O(n)$。

為了加速運算，目標是要找到 $\mod p$ 下的循環長度 $L$，最後求出$F_{2^n \mod L}$，根據待會的證明，保證 $L \le p$，那複雜度就可以回到 $O(\log L)$ 解決。但為了要找到 $L$ 又是一段很長的故事，總時間複雜度為 $O(\sqrt{p})$，不用保證 $p$ 是質數。

參考資料

故事

數列$F_0 = 1, F_1 = 1, F_2 = 2, \cdots$，循環是連續兩項出現重複，而費氏數列會完全循環，也就是出現連續兩項$F_i = 0, F_{i-1} = 1$。下方是一個 $\mod 4$ 的情況。

1
1 1 2 3 1 0 | 1 2 3 1 0 | ...

要找到恰好連續兩項$F_i = 0, F_{i-1} = 1$ 是困難的，考慮去找到$F_i = 0$ 即可，接著再去想辦法讓$F_{i-1} = 1$。

假設最小的 $k$ 滿足$F_k = 0 \mod p$，而$F_{k-1} = a \mod p$，那麼之後的序列$F_{i} = a^j F_{i+j \times k} \mod p$。從矩陣乘法的概念中可以理解，是一個常數為 $a$ 的初始項，第二輪循環常數就會變成 $a^2$，類推。

接下來

  • 考慮一個嚴重的問題「 何種模 $p$ 情況一定循環，即從$F_0$ 再次循環。 」答案是 質數 $p$ 。
    原因是$F_{i} = a^j F_{i+j \times k} \mod p$，由於 $a$ 與 $p$ 互質，$a^j \mod p \neq 0$ 恆成立，同時還是一個 $ord_{p}(a) = p-1$，這部分從歐拉定理中可以了解，那麼只有可能在$F_{i} = 0$ 的情況成立，就是 $k$ 的倍數之外，不發生$F_i = 0$ 的出現。
  • 接續上一個問題「模 $p$ 不是質數怎麼處理？」
    進行質因數分解，對於每一個質因子找到模循環長度，模 $p$ 循環長度就是所有質因子循環長度的最小公倍數 lcm。

現在問題落在 $k$ 怎麼找到，若能找到 $k$，其循環長度落在 $k$ 的倍數，或者有更好的獲取方式。

  • 若模質數 $p$ 且滿足 $p > 5$，5 是 $p$ 的二次剩餘 (quadratic residue)，意即滿足 $\exists \; x^2 \equiv 5 \mod p$，循環長度為 $p-1$ 的因數。反之，循環長度是 $2(p+1)$ 的因數。

關於二次剩餘的判斷，在模質數 $p$ 下，對於 $gcd(x, p) = 1$，藉由歐拉定理得到 $x^{p-1} \equiv 1 \mod p$，以下不保證是正確的說法，提供理解的一個方案。

  • 若 $d$ 是 $p$ 的二次剩餘，則滿足 $d^{(p-1)/2} \equiv 1 \mod p$，因為$x^{2 \times (p-1)/2} \equiv 1 \mod p \Rightarrow x^{p-1} \equiv 1 \mod p$。
  • 若非二次剩餘，則滿足 $d^{(p-1)/2} \equiv -1 \mod p$，因為 $\left [ d^{(p-1)/2} \right ]^2 \equiv 1 \mod p \Rightarrow d^{p-1} \equiv 1 \mod p$
  • 數學上用 Legendre symbol 來表示這個判斷 wiki。

回過頭來，看一下費氏數列的公式解

$F_n = \frac{1}{\sqrt{5}} \left [ \left ( \frac{1+\sqrt{5}}{2} \right )^n - \left (\frac{1-\sqrt{5}}{2} \right )^n \right ]$

藉由展開公式，儘管它是實數、根號，展開之後一定只會剩下整數冪次的總和。令 $a = \sqrt{5}$，觀察二次剩餘與否和滿足$F_n = 0$ 的關係。

在模質數 $p$ 下，滿足二次剩餘$F_n \equiv 0 \mod p$，當 $n = p-1$ 的時候成立，可以藉由噁心的展開式得到。同理在非二次剩餘情況，$n = 2(p+1)$，找到一個最大的倍數情況，答案一定落在其因數下。詳細推導請看參考資料，太噁心就不提。

參考資料中有特別提到，有一個地方還 沒有確認 ，對於模數 $p^k$ 的循環長度 $g(p) \times p^{k-1}$ 如何證明。但我想根據中國餘式定理，能了解循環長度倍數的模關係吧。接著由於大整數分解期望是 $O(\sqrt{n})$，中間也要找到所有因數來得到循環長度的驗證，還要搭配快速矩陣乘法，最後也是 $O(\sqrt{n})$。

這一題是更進階的費氏數列計算，從建表、矩陣乘法，最後來到循環搜尋，這些證明不久之後就會忘記了，二次剩餘判斷的想法還會留著，其他就忘得一乾二淨吧。

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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
#include <bits/stdc++.h> 
using namespace std;
#define MILLER_BABIN 4
typedef unsigned long long UINT64;
UINT64 mul(UINT64 a, UINT64 b, UINT64 mod) { 
    UINT64 ret = 0; 
    for (a = a >= mod ? a%mod : a, b = b >= mod ? b%mod : b; b != 0; b>>=1, a <<= 1, a = a >= mod ? a - mod : a) { 
        if (b&1) {
            ret += a;
            if (ret >= mod) 
                ret -= mod;
        } 
    } 
    return ret; 
}
struct Matrix {
    UINT64 v[2][2];
    int row, col; // row x col
    Matrix(int n, int m, int a = 0) {
        memset(v, 0, sizeof(v));
        row = n, col = m;
        for(int i = 0; i < row && i < col; i++)
            v[i][i] = a;
    }    
    Matrix multiply(const Matrix& x, const long long mod) const {
        Matrix ret(row, x.col);
        for(int i = 0; i < row; i++) {
            for(int k = 0; k < col; k++) {
                if (!v[i][k])
                	continue;
                for(int j = 0; j < x.col; j++) {
                    ret.v[i][j] += mul(v[i][k], x.v[k][j], mod);
                    if (ret.v[i][j] >= mod)
                    	ret.v[i][j] -= mod;
                }
            }
        }
        return ret;
    }
    Matrix pow(const long long& n, const long long mod) const {
        Matrix ret(row, col, 1), x = *this;
        long long y = n;
        while(y) {
            if(y&1)	ret = ret.multiply(x, mod);
            y = y>>1, x = x.multiply(x, mod);
        }
        return ret;
    }
} FibA(2, 2, 0);
#define MAXL (50000>>5)+1
#define GET(x) (mark[x>>5]>>(x&31)&1)
#define SET(x) (mark[x>>5] |= 1<<(x&31))
int mark[MAXL], P[50000], Pt = 0;
void sieve() {
    register int i, j, k;
    SET(1);
    int n = 46340;
    for (i = 2; i <= n; i++) {
        if (!GET(i)) {
            for (k = n/i, j = i*k; k >= i; k--, j -= i)
                SET(j);
            P[Pt++] = i;
        }
    }
}
UINT64 mpow(UINT64 x, UINT64 y, UINT64 mod) { // mod < 2^32 
    UINT64 ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
UINT64 mpow2(UINT64 x, UINT64 y, UINT64 mod) {
    UINT64 ret = 1;
    while (y) {
        if (y&1)
            ret = mul(ret, x, mod);
        y >>= 1, x = mul(x, x, mod);
    }
    return ret;
} 
void exgcd(long long x, long long y, long long &g, long long &a, long long &b) {
    if (y == 0)
        g = x, a = 1, b = 0;
    else
        exgcd(y, x%y, g, b, a), b -= (x/y) * a;
}
long long llgcd(long long x, long long y) {
    if (x < 0)    x = -x;
    if (y < 0)    y = -y;
    if (!x || !y)    return x + y;
    long long t;
    while (x%y)
        t = x, x = y, y = t%y;
    return y;
}
long long inverse(long long x, long long p) {
    long long g, b, r;
    exgcd(x, p, g, r, b);
    if (g < 0)	r = -r;
    return (r%p + p)%p;
}
int isPrime(long long p) { // implements by miller-babin
    if (p < 2 || !(p&1))	return 0;
    if (p == 2)				return 1;
    long long q = p-1, a, t;
    int k = 0, b = 0;
    while (!(q&1))	q >>= 1, k++;
    for (int it = 0; it < MILLER_BABIN; it++) {
        a = rand()%(p-4) + 2;
        t = mpow2(a, q, p);
        b = (t == 1) || (t == p-1);
        for (int i = 1; i < k && !b; i++) {
            t = mul(t, t, p);
            if (t == p-1)
                b = 1;
        }
        if (b == 0)
            return 0;
    }
    return 1;
}
long long pollard_rho(long long n, long long c) {
    long long x = 2, y = 2, i = 1, k = 2, d;
    while (true) {
        x = (mul(x, x, n) + c);
        if (x >= n)	x -= n;
        d = llgcd(x - y, n);
        if (d > 1) return d;
        if (++i == k) y = x, k <<= 1;
    }
    return n;
}
void factorize(int n, vector<long long> &f) {
    for (int i = 0; i < Pt && P[i]*P[i] <= n; i++) {
    	if (n%P[i] == 0) {
    		while (n%P[i] == 0)
    			f.push_back(P[i]), n /= P[i];
    	}
    }
    if (n != 1)	f.push_back(n);
}
void llfactorize(long long n, vector<long long> &f) {
    if (n == 1)	
        return ;	
    if (n < 1e+9) {
        factorize(n, f);
        return ;
    }
    if (isPrime(n)) {
        f.push_back(n);
        return ;
    }
    long long d = n;
    for (int i = 2; d == n; i++)
        d = pollard_rho(n, i);
    llfactorize(d, f);
    llfactorize(n/d, f);
}
// above largest factor
// ---------------------- //
int legendre_symbol(UINT64 d, UINT64 p) {
    if (d%p == 0)	return 0;
    return mpow2(d, (p-1)>>1, p) == 1 ? 1 : -1;
}
void factor_gen(int idx, long long x, vector< pair<long long, int> > &f, vector<long long> &ret) {
    if (idx == f.size()) {
        ret.push_back(x);
        return ;
    }
    for (long long i = 0, a = 1; i <= f[idx].second; i++, a *= f[idx].first)
        factor_gen(idx+1, x*a, f, ret);
}
void factor_gen(long long n, vector<long long> &ret) {
    vector<long long> f;
    vector< pair<long long, int> > f2;
    llfactorize(n, f);
    sort(f.begin(), f.end());
    int cnt = 1;
    for (int i = 1; i <= f.size(); i++) {
        if (i == f.size() || f[i] != f[i-1])
            f2.push_back(make_pair(f[i-1], cnt)), cnt = 1;
        else
            cnt ++;
    }
    factor_gen(0, 1, f2, ret);
    sort(ret.begin(), ret.end());
}
UINT64 cycleInFib(UINT64 p) {
    if (p == 2)	return 3;
    if (p == 3)	return 8;
    if (p == 5)	return 20;
    vector<long long> f;
    if (legendre_symbol(5, p) == 1)
        factor_gen(p-1, f);
    else
        factor_gen(2*(p+1), f);
    long long f1, f2;
    for (int i = 0; i < f.size(); i++) {
        Matrix t = FibA.pow(f[i]-1, p);
        f1 = (t.v[0][0] + t.v[0][1])%p;
        f2 = (t.v[1][0] + t.v[1][1])%p;
        if (f1 == 1 && f2 == 0)
            return f[i];
    }
    return 0;
}
UINT64 cycleInFib(UINT64 p, int k) {
    UINT64 s = cycleInFib(p);
    for (int i = 1; i < k; i++)
        s = s * p;
    return s;
}
int main() {
    sieve();
    FibA.v[0][0] = 1, FibA.v[0][1] = 1;
    FibA.v[1][0] = 1, FibA.v[1][1] = 0;
    
    int testcase;
    scanf("%d", &testcase);
    while (testcase--) {
        long long n, m;
        scanf("%lld %lld", &n, &m);
        
        vector<long long> f;
        map<long long, int> r;
            
        llfactorize(m, f);
        for (auto &x : f)
            r[x]++;
            
        UINT64 cycle = 1; 
        for (auto &x : r) {
            UINT64 t = cycleInFib(x.first, x.second);
            cycle = cycle / llgcd(t, cycle) * t;
        }
        
        n = mpow2(2, n, cycle);
        Matrix t = FibA.pow(n, m);
        long long fn = t.v[1][0];
        printf("%lld\n", fn);
    }
    return 0;
}