原题链接:P4000 斐波那契数列
(一开始以为是道水题,没想到竟然如此复杂,呜呜呜)
题目:求 f_n\bmod P 的值。
注: 原题的模数是小写的 p,此处为了区分下文的质数 p,改为大写的 P。
首先,n 达到了 10^{30000000} 的级别,矩阵快速幂肯定 TLE。
于是我们想到了要通过一些方法降低 n 的大小,最先想到的是对 n 模以某个数,然后再矩阵快速幂不就好了嘛。
那么先手算一下取模的这个数和 P 的关系吧,假设 P=4,那么斐波那契数列长这样:
1,1,2,3,1,0,1,1,2,3,1,0,...
好像确实有规律,针对 P=4 时我们只需求 f_{n\bmod 6} 就好了!
可是对于其他的 P 要怎么求?为了解决这个问题,我们可能需要一点前置知识……
前置知识:二次剩余
可能需先掌握「原根」的知识。
二次剩余是指这样一类问题:
已知 n,p 的值,解方程 x^2\equiv n\pmod p 。
洛谷上有模板题:P5491 【模板】二次剩余。此处不做详述,只讲与下文有关的,即 p 是奇素数的情况。
首先 p\mid n 时,x^2\equiv n\equiv 0\pmod p,故 x=0。
当 p\nmid n 时,x 可能无解,如 n=2,p=3 时。那么 n,p 满足什么条件时 x 有解呢?
结论一:n^{\frac{p-1}{2}}\equiv 1\pmod p 是 x 有不为 0 的解的充要条件。以下证明:
当 x 有不为 0 的解时,x^{p-1}\equiv 1。因此 n^{\frac{p-1}{2}}\equiv(x^2)^{\frac{p-1}{2}}\equiv x^{p-1}\equiv 1,必要性证毕。
当 n^{\frac{p-1}{2}}\equiv 1 时,设 g 是 p 的原根,则存在 k 使得 n=g^k,故 g^{\frac{k}{2}(p-1)}\equiv 1。因为 g 的阶是 p-1,所以 (p-1)|\frac{k}{2}(p-1),即 k 是偶数。那么解得 x\equiv g^{\frac{k}{2}},有不为 0 的解。充分性证毕。
那什么时候无解呢?因为 n^{p-1}\equiv 1,(n^{\frac{p-1}{2}}-1)(n^{\frac{p-1}{2}}+1)\equiv 0,即 n^{\frac{p-1}{2}}\equiv \pm1。由结论一推理得,n^{\frac{p-1}{2}}\equiv -1 是 x 无解的充要条件,记为结论二。
事实上,上面这些对下文有用的只有这句话:
------------
## 斐波那契循环节
------------
接下来就是主线情节了。
**结论三**:斐波那契数列 $\{f_i\}$ 对 $p$ 取模后必定是周期数列。
显然,只要存在 $i\neq j$ 满足 $f_i\equiv f_j$ 且 $f_{i+1}\equiv f_{j+1}$。那么由递推关系可以得到 $f_{i+k}\equiv f_{j+k}$,命题就成立了。
要证明存在这样的 $i,j$ 也不难,因为不可能有无穷多对不同的 $(f_i,f_{i+1})$,至多 $p^2$ 对,因此一定存在 $i\neq j$,满足 $f_i\equiv f_j$ 且 $f_{i+1}\equiv f_{j+1}$。故循环节一定存在,且长度不超过 $p^2$,证毕。
设 $\pi(p)$ 是这个数列的周期(注意,此处并不要求 $\pi(p)$ 是最小正周期)。那要怎样求 $\pi(p)$?
不难发现,若存在 $k$ 满足 $\forall i, f_{i}\equiv f_{i+k}\pmod p$,那 $k$ 就是 $\{f_i\bmod p\}$ 的周期,取 $\pi(p)=k$ 就好了。
------------
### Part 1:
让我们从最简单的情况开始:当模数 $P$ 等于**某个质数** $p$ 的时候。
先搬出斐波那契数列的通项公式:$f_i=\dfrac{(\frac{1+\sqrt 5}{2})^i-(\frac{1-\sqrt 5}{2})^i}{\sqrt 5}$。
因为公式里有 $\sqrt 5,\frac1 2$,所以我们先单独讨论 $p=2,5$ 的情况:经计算取 $\pi(2)=3,\pi(5)=20$。
撇除这些特殊情况后,$p$ 便是一个不等于 $5$ 的奇素数。但根据二次剩余的知识,$\sqrt 5$ 模 $p$ 时不一定能表示成整数,需分类讨论。
- $x^2\equiv 5\pmod p$ 有非 $0$ 解
这时 $x$ 等价于 $\sqrt 5$ ,那么 $\sqrt 5,\frac{1}{\sqrt 5},\frac 1 2$ 模 $p$ 就都有意义,即 $\frac{1+\sqrt 5}{2},\frac{1-\sqrt 5}{2}$ 能用整数表示,适用于欧拉定理。所以:
$$\begin{aligned}f_{i+\varphi(p)}&\equiv \frac{(\frac{1+\sqrt 5}{2})^{i+\varphi(p)}-(\frac{1-\sqrt 5}{2})^{i+\varphi(p)}}{\sqrt 5}\\&\equiv\frac{(\frac{1+\sqrt 5}{2})^i-(\frac{1-\sqrt 5}{2})^i}{\sqrt 5}=f_i\end{aligned}$$
故 $\varphi(p)$ 是数列的一个周期,取 $\pi(p)= \varphi(p)=p-1$。
- $x^2\equiv 5\pmod p$ 无解
无解意味着找不到整数替代 $\sqrt 5$,需要**扩域**。
扩域指的是令某个数 $w$ 满足 $w^2\equiv K$,其中 $K^{\frac{p-1}{2}}\equiv -1$。因为 $x^2\equiv K$ 无解,所以这里的 $w$ 类似于 $i^2=-1$ 的 $i$,并非整数。同时,由 $w$ 构成的数 $a+bw$ 也不满足欧拉定理等,需要新的结论。
**结论四**:对于奇素数 $p$,有 $(a+bw)^{p+1}\equiv a^2-b^2K\pmod p$。
以下证明:由二项式定理可得:
$$(a+bw)^{p+1}=\sum_{i=0}^{p+1}C_{p+1}^ia^{p+1-i}(bw)^i$$
聚焦于组合数 $C_{p+1}^i\equiv\frac{p!}{i!(p+1-i)!}$,可以发现仅在 $i=0,1,p,p+1$ 这四个值时组合数模 $p$ 后不为 $0$。因此原式等于:
$$C_{p+1}^0a^{p+1-0}(bw)^0+C_{p+1}^1a^{p+1-1}(bw)^1+C_{p+1}^{p}a^{p+1-p}(bw)^p+C_{p+1}^{p+1}a^{p+1-p-1}(bw)^{p+1}$$
即 $a^{p+1}+a^{p}bw+ab^pw^p+b^{p+1}w^{p+1}\equiv a^2+abw+abw^p+b^2w^{p+1}$。
$\because w^2\equiv K,K^{\frac{p-1}{2}}\equiv -1
\therefore (w^2)^{\frac{p-1}{2}}\equiv w^{p-1}\equiv -1
代入原式得:a^2-b^2w^2\equiv a^2-b^2K,证毕。
接下来,设 w^2\equiv 5,则通项公式可写成:f_i=\dfrac{(\frac{p+1}{2}+\frac{p+1}{2}w)^i-(\frac{p+1}{2}-\frac{p+1}{2}w)^i}{w}。
先看 \frac{p+1}{2}+\frac{p+1}{2}w,由结论四得:(\frac{p+1}{2}+\frac{p+1}{2}w)^{p+1}\equiv (\frac{p+1}{2})^2-5(\frac{p+1}{2})^2\equiv -1。
所以 (\frac{p+1}{2}+\frac{p+1}{2}w)^{2p+2}\equiv 1,同理 (\frac{p+1}{2}-\frac{p+1}{2}w)^{2p+2}\equiv 1。
因此 f_{i+2p+2}\equiv f_i,取 \pi(p)=2p+2。
Part 2:
接下来,讨论模数 P 是某个质数的幂的情况,记 P=p^k。
注:Part 2 的证明参考了部分资料(链接见文末)。
结论五:对于质数 p,已知 a\equiv 1\pmod p,则 a^{p^k}\equiv 1\pmod {p^{k+1}}。
以下证明:考虑数学归纳法,假设 k 时等式成立,证明 k+1 时等式也成立。
设 a^{p^{k}}=sp^{k+1}+1,那么 a^{p^{k+1}}=(a^{p^k})^p=(sp^{k+1}+1)^p。
由二项式定理:
(sp^{k+1}+1)^p=\sum_{i=0}^pC_p^i(sp^{k+1})^i
当 i\ge 2 时,p^{(k+1)i}\bmod p^{k+2}=0,只需算出 i=0,1 的值。代入得 1+sp^{k+2}\equiv 1\pmod {p^{k+2}},证毕。
方便起见,记 \frac{1+\sqrt 5}{2}=u,\frac{1-\sqrt 5}{2}=v。
结论六:对于质数 p,有 u^{\pi(p)}\equiv v^{\pi(p)}\equiv 1\pmod p。
\because f_{\pi(p)}\equiv 0\pmod p
\therefore u^{\pi(p)}\equiv v^{\pi(p)}
\because f_{\pi(p)+1}\equiv 1
\therefore {u^{\pi(p)}u-v^{\pi(p)}v}\equiv \sqrt 5
代入 u^{\pi(p)}\equiv v^{\pi(p)} 化简得:u^{\pi(p)}\equiv v^{\pi(p)}\equiv 1,证毕。
结合结论五和结论六可得:u^{\pi(p)p^{k-1}}\equiv v^{\pi(p)p^{k-1}}\equiv 1\pmod {p^k}。
因此 f_{\pi(p)p^{k-1}}\equiv \frac{1}{\sqrt 5}(u^{\pi(p)p^{k-1}}-v^{\pi(p)p^{k-1}})\equiv 0\pmod {p^k}。
且 f_{1+\pi(p)p^{k-1}}\equiv \frac{1}{\sqrt 5}(u^{\pi(p)p^{k-1}}u-v^{\pi(p)p^{k-1}}v)\equiv \frac{u-v}{\sqrt 5}\equiv 1\pmod {p^k}。
也就是说:f_{i+\pi(p)p^{k-1}}\equiv f_i\pmod {p^k},取 \pi(p^k)=\pi(p)p^{k-1}。
Part 3:
最后一步了,也是最简单的情况。
结论七:设模数 P 的质因数分解式为 \prod p^k,则可取 \pi(P)=\operatorname{lcm}(\pi(p_i^{k_i}))。
以下证明:
\because f_{\pi(P)}\equiv 0\pmod{p_i^{k_i}},p_i^{k_i}\mid f_{\pi(P)}
\therefore {\pi(p_i^{k_i})}\mid {\pi(P)}
取 \pi(P) 为 \pi(p_i^{k_i}) 的最小公倍数即可,证毕。
代码实现
回到原题,开始整理思路。
-
对 P 进行质因数分解,时间复杂度 O(\sqrt P);
-
对于 P 的每个质因数分解 p_i^{k_i},求出对应的周期长 \pi(p_i^{k_i});
-
求出每个周期的最小公倍数,作为 \pi(P) 的值;
-
把读入的 n 模以 \pi(P),然后矩阵快速幂解决。
就这样,代码在此:
#include <iostream>
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;
ll FastPow(ll a, ll b, ll mod){
ll res = 1;
while(b){
if(b & 1) res = res * a % mod;
b >>= 1, a = a * a % mod;
}
return res;
}
ll FastPowNoMod(ll a, ll b){
ll res = 1;
while(b){
if(b & 1) res = res * a;
b >>= 1, a = a * a;
}
return res;
}
const int Mx = 3e7 + 10;
int Fac[30], Tim[30], tot;
char N[Mx];
void Divide(int p){ // 分解质因数
for(int i = 2; i * i <= p; ++i) if(p % i == 0){
Fac[++tot] = i;
while(p % i == 0) p /= i, ++Tim[tot];
}
if(p != 1) Fac[++tot] = p, Tim[tot] = 1;
}
ll PrimeLoop(ll p){ // 求 pi(p)
if(p == 2) return 3;
if(p == 5) return 20;
if(FastPow(5, (p - 1) >> 1, p) == 1) return p - 1;
return 2 * p + 2;
}
ll PrimePow(int id){ return FastPowNoMod(Fac[id], Tim[id] - 1) * PrimeLoop(Fac[id]);} // 求 pi(p^k)
ll Gcd(ll a, ll b){ return b == 0 ? a : Gcd(b, a % b);}
ll GetLoop(){ // 求 pi(P)
ll res = PrimePow(1), temp;
for(int i = 2; i <= tot; ++i) temp = PrimePow(i), res = res / Gcd(res, temp) * temp;
return res;
}
int P;
struct Matrix{
ll val[3][3];
Matrix operator * (const Matrix& a) const {
Matrix Res = (*this);
for(int i = 1; i <= 2; ++i) for(int j = 1; j <= 2; ++j)
Res.val[i][j] = (val[i][1] * a.val[1][j] % P + val[i][2] * a.val[2][j] % P) % P;
return Res;
}
Matrix operator *= (const Matrix& a) { return (*this) = a * (*this);}
Matrix operator ^ (ll k) const {
Matrix Res, Temp = (*this);
Res.val[1][1] = Res.val[2][2] = 1, Res.val[1][2] = Res.val[2][1] = 0;
while(k){
if(k & 1) Res *= Temp;
k >>= 1, Temp *= Temp;
}
return Res;
}
} Fibo, Cell;
void Solve(ll n){
if(n == 0){
printf("0");
return;
}
Fibo.val[1][1] = 0, Fibo.val[2][1] = Fibo.val[1][2] = Fibo.val[2][2] = 1;
Cell.val[1][1] = 1, Cell.val[1][2] = 1;
Matrix Res = Cell * (Fibo ^ (n - 1));
printf("%lld", Res.val[1][1] % P);
}
int main(){
scanf(" %s %d", N, &P);
if(P == 1){
printf("0");
return 0;
}
Divide(P);
ll Loop = GetLoop(), n = 0;
int Len = strlen(N);
for(int i = 0; i < Len; ++i) n = (n * 10 + (N[i] ^ 48)) % Loop;
Solve(n);
return 0;
}
参考资料: