Description
农夫栋栋近年收入不景气, 正在他发愁如何能多赚点钱时, 他听到隔壁的小朋友在讨论兔子 繁殖的问题.
问题是这样的: 第一个月初有一对刚出生的小兔子, 经过两个月长大后, 这对兔子从第三个 月开始, 每个月初生一对小兔子. 新出生的小兔子生长两个月后又能每个月生出一对小兔子. 问第 \(n\) 个月有多少只兔子?
聪明的你可能已经发现, 第 \(n\) 个月的兔子数正好是第 \(n\) 个 Fibonacci (斐波那契) 数. 栋栋不懂什么是 Fibonacci 数, 但他也发现了规律: 第 \(i+2\) 个月的兔子数等于第 \(i\) 个月 的兔子数加上第 \(i+1\) 个月的兔子数. 前几个月的兔子数依次为:
\[1, 1, 2, 3, 5, 8, 13, 21, 34, \ldots\]
栋栋发现越到后面兔子数增长的越快, 期待养兔子一定能赚大钱, 于是栋栋在第一个月初 买了一对小兔子开始饲养.
每天, 栋栋都要给兔子们喂食, 兔子们吃食时非常特别, 总是每 \(k\) 对兔子围成一圈, 最后 剩下的不足 \(k\) 对的围成一圈, 由于兔子特别害怕孤独, 从第三个月开始, 如果吃食时围成某一个圈的只有一对兔子, 这对兔子就会很快死掉.
我们假设死去的总是刚出生的兔子, 那么每个月的兔子数仍然是可以计算的. 例如, 当 \(k = 7\) 时, 前几个月的兔子数依次为:
\[1, 1, 2, 3, 5, 7, 12, 19, 31, 49, 80, \ldots\]
给定 \(n\), 你能帮助栋栋计算第 \(n\) 个月他有多少对兔子么? 由于答案可能非常大, 你只 需要告诉栋栋第 \(n\) 个月的兔子对数除 \(p\) 的余数即可.
Input
输入一行, 包含三个正整数 \(n, k, p\).
Output
输出一行, 包含一个整数, 表示栋栋第 \(n\) 个月的兔子对数除 \(p\) 的余数.
Sample Input
6 7 100Sample Output
7Data Range
对于所有数据, 满足 \(1 \leq n \leq 10^{18}, 2 \leq k \leq 10^6, 2 \leq p \leq 10^9\)
Explanation
其实 vfk 说的很详细了, 原文链接:http://vfleaking.blog.163.com/blog/static/174807634201341721051604/
首先我们考察模 \(G\) 意义下的 Fibonacci 数列的循环节吧～
这个有个高端的名字叫 Pisano period.
可以去看 http://math.arizona.edu/~ura-reports/071/Campbell.Charles/Final.pdf (下载链接:pisano_period.pdf
超好懂.
OEIS 列出了前几项:http://oeis.org/A001175
现在我们来证, 循环节 \(p \leq 6G\)
设 \(G = 2^t \cdot 5^s \cdot p_1^{k_1} \cdot p_2^{k_2} \cdot \ldots \cdot p_a^{k_a}\)
\(t = 0\) 或 \(s = 0\) 的情况就不考虑了.
\[P \leq lcm(3 \cdot 2^{t - 1}, 4 \cdot 5^s, 2(p_1 + 1) \cdot p_1^{k_1 - 1}, 2(p_2 + 1) \cdot p_2^{k_2 - 1}, \ldots, 2(p_a + 1) \cdot p_a^{k_a - 1})\]
\[= lcm(3, 2^{t - 1}, 4, 5^s, \frac{p_1 + 1}{2} \cdot p_1^{k_1 - 1}, \frac{p_2 + 1}{2} \cdot p_2^{k_2 - 1}, \ldots, \frac{p_a + 1}{2} \cdot p_a^{k_a - 1})\]
\[\leq 3 \cdot 4 \cdot (2^{t - 1} \cdot 5^s \cdot p_1^{k_1} \cdot p_2^{k_2} \cdot \ldots \cdot p_a^{k_a})\]
\[\leq 12 \cdot \frac{G}{2}\]
\[= 6G\]
这个是不可改进的, 因为形如 \(2 \cdot 5^s\) 的数能达到这个上界.
……好帅气的结论!! !
记 Fibonacci 数列第 \(n\) 项对 \(G\) 取模的值为 \(Fib[n]\).
考虑题目中给的 \(F\) 与 \(Fib\) 的区别, 就在于那个减一了～～
神马时候会减一呢?
考虑第一次减一, 假设出现在 \(k\):
则 \(Fib[k] = 1\)
那么此时:
...... Fib [k-2] Fib [k-1] 0 Fib [k-1] Fib [k-1] 2Fib [k-1] 3Fib [k-1]. ......
考虑后面那一串, 现在暂时不对 G 取模鸟～:
Fib [k-1] Fib [k-1] 2Fib [k-1] 3Fib [k-1] 5Fib [k-1] 8Fib [k-1]. ......
它们都是 \(Fib[k - 1]\) 的倍数耶……
同时除以 \(Fib[k - 1]\) 就有:
1 1 2 3 5 8......
新的 Fibonacci 数列.
于是我们注意到每次减一之后一定会出现一个 0 咯……而 0 前面一个数如果是 a, 那么这一段的数列就是:
...... a 0 Fib [1] a Fib [2] a Fib [3] a Fib [4] a......
我们把 “Fib [1] a Fib [2] a Fib [3] a Fib [4] a...... Fib [L-1] a 0” 这样的一串称作数列 b [a].
那么数列 b 的长度是多少呢? 就是问 L 等于多少, 就是问使得 Fib [L] a=1 (mod G) 的最小的 L.
由于 Fib 循环节最长是 \(6G \leq 6 \cdot 10^6\), 所以我们可以暴力记录下来 Fib 的一个循环节, 用逆元乱搞下就能求出 L.
我们还能求出 \(next[a]\) 数组, 表示 “\(Fib[L - 1]a\)” 的值.
这样如果 \(b[a]\) 数列出现了, 后面一定跟着 \(b[next[a]]\) 数列.
当然如果不存在 \(a\) 的逆元或者没有 \(Fib[L] = a\) 的逆元, 就永远不会有下一个 \(b\) ……
此时反正不会减一, 直接当 Fibonacci 数列处理.
用矩阵乘法算模 \(P\) 的答案.
暴力搞搞都行了.
Source Code
#include <iostream>
#include <cstdlib>
#include <cstdio>
#define rep(_var,_begin,_end) for(int _var=_begin;_var<=_end;_var++)
using namespace std;
typedef long long lli;
const int maxn = 1001000;
const lli infinit = 0x7fffffffffffffffll;
int modr = 1000000007;
struct matrix { lli arr[3][3];
lli& operator() (int x, int y) { return this->arr[x][y]; }
lli operator() (int x, int y) const { return this->arr[x][y]; }
matrix& operator *= (const matrix &b) {
lli c[3][3];
rep(i, 0, 2) rep(j, 0, 2) {
c[i][j] = 0;
rep(k, 0, 2) c[i][j] += arr[i][k] * b.arr[k][j];
c[i][j] %= modr; }
rep(i, 0, 2) rep(j, 0, 2)
this->arr[i][j] = c[i][j];
return *this;
}
matrix(lli _1, lli _2, lli _3, lli _4, lli _5, lli _6, lli _7, lli _8, lli _9) {
arr[0][0] = _1, arr[0][1] = _2, arr[0][2] = _3;
arr[1][0] = _4, arr[1][1] = _5, arr[1][2] = _6;
arr[2][0] = _7, arr[2][1] = _8, arr[2][2] = _9;
return ;
}
matrix(void) {
rep(i, 0, 2) rep(j, 0, 2)
arr[i][j] = 0;
return ;
}
matrix(int k) {
rep(i, 0, 2) rep(j, 0, 2)
arr[i][j] = i == j ? k : 0;
return ;
}
void print(void) {
rep(i, 0, 2) { rep(j, 0, 2)
printf("%lld ", arr[i][j]);
printf("\n"); }
return ;
}
};
matrix operator ^ (matrix a, lli pw) {
matrix res(1);
for (lli i = pw; i > 0; i >>= 1) {
if (i & 1) res *= a;
a *= a;
}
return res;
}
lli m_rev(lli a, lli m)
{
lli x1 = 1, x2 = 0, x3 = a;
lli y1 = 0, y2 = 1, y3 = m;
// Extended euclid algorithm.
while (y3 != 0) {
lli k = x3 / y3;
lli t1 = x1 - y1 * k, t2 = x2 - y2 * k, t3 = x3 - y3 * k;
x1 = y1, x2 = y2, x3 = y3;
y1 = t1, y2 = t2, y3 = t3;
}
if (x3 != 1)
return 0;
return x1 >= 0 ? x1 : x1 + m;
}
lli n, K, p;
int first_pos[maxn], fib[maxn*6+1];
lli rev[maxn], pow[maxn];
int next[maxn];
bool vis[maxn + 1];
int main(int argc, char** argv)
{
scanf("%lld%lld%lld", &n, &K, &p);
modr = p;
// Calculate repetition subsequence
for (int i = 0; i < maxn; i++)
first_pos[i] = -1;
fib[0] = fib[1] = 1;
for (int i = 2; true; i++) {
int cur = fib[i-1] + fib[i-2];
while (cur > K)
cur -= K;
fib[i] = cur;
if (first_pos[cur] < 0)
first_pos[cur] = i;
if (fib[i] == 1 && fib[i-1] == 1)
break;
}
// Multiplicative inverse calculation
for (int i = 0; i < K; i++)
rev[i] = m_rev(i, K);
for (int i = 0; i < K; i++) {
if (!rev[i] || !first_pos[rev[i]]) {
next[i] = -1, pow[i] = infinit;
} else {
int k = first_pos[rev[i]];
next[i] = (lli)fib[k-1] * i % K;
pow[i] = k+1;
}
}
// Matrix multiplication
matrix matA(1, 1, 0,
1, 0, 0,
0, 0, 1);
matrix matB(1, 1, 0,
1, 0, 0,
modr - 1, 0, 1);
int cirN;
for (cirN = 1; cirN != -1 && !vis[cirN]; cirN = next[cirN])
vis[cirN] = true;
matrix matRes(1);
for (int i = 1; i != cirN && n > 0; i = next[i]) {
if (n < pow[i]) {
matRes *= matA ^ n;
n = 0;
} else {
matRes *= matA ^ (pow[i] - 1);
matRes *= matB;
n -= pow[i];
}
}
// Final processing
if (n > 0) {
int cur = cirN;
lli tot = 0;
matrix matComp(1);
do {
matComp *= matA ^ (pow[cur] - 1);
matComp *= matB;
tot += pow[cur];
cur = next[cur];
} while (cur != cirN);
matRes *= matComp ^ (n / tot);
n %= tot;
for (int i = cirN; n > 0; i = next[i]) {
if (n < pow[i]) {
matRes *= matA ^ n;
n = 0;
} else {
matRes *= matA ^ (pow[i] - 1);
matRes *= matB;
n -= pow[i];
}
}
}
lli vec[3] = {0, 1, 1};
lli res = 0;
rep(k, 0, 2) res += vec[k] * matRes(k, 0);
res %= modr;
printf("%lld\n", res);
return 0;
}