我居然没看题解瞎搞出来了？
题解：
不难想到找到每次减1的位置，然后减去它对最终答案的贡献。
假设有一个地方是(x,1(mod~k))
那么减了1后就变成了(x,0)。
然后可以推到(x,0,x,x)。
可以看做以(x,x)为开头，做新的序列。
设y为x在mod k意义下的逆元，那么下次1的地方就是原斐波拉契序列中第一次出现y的位置。
如果没有逆元那就结束了。
原斐波拉契序列的非循环长度是(O(k))级别的，这个可以通过随机序列第一次出现相同元素来理解，为在mo意义下应该是可以看做随机的，值域大小是(O(k^2))，那么期望长度就是(O(sqrt {k^2}))。
所以斐波拉契序列做个3k左右然后预处理每一个元素第一次出现的位置。
然后x也是有循环的，因为x<k，所以肯定是(O(k))的，那么就可以一起做了。
最后就是个等比矩阵求和，由于不一定有逆，所以要用分治法去求等比矩阵和。
由于循环开始之前和最后结尾的地方都要判，所以有点复杂。
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i < B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
ll n; int k, mo;
const int N = 3e6 + 5;
int f[N], fi[N];
int p[N], p0, us[N];
int gcd(int x, int y) { return (!y ? x : gcd(y, x % y));}
void exgcd(int a, int b, int &x, int &y) {
if(!b) {x = a, y = 0; return;}
exgcd(b, a % b, y, x); y -= (a / b) * x;
}
int qni(int p, int q) {
int x, y;
exgcd(p, q, x, y);
x = (x % q + q) % q;
return x;
}
int l = -1, r;
void work() {
int x = 1;
while(1) {
if(gcd(x, k) > 1) break;
int y = qni(x, k);
if(fi[y]) {
p0 ++;
p[p0] = p[p0 - 1] + fi[y];
if(us[y]) {
l = us[y]; r = p0;
break;
}
us[y] = p0;
x = (ll) f[fi[y] - 1] * x % k;
} else break;
}
}
struct jz {
ll a[2][2];
jz() {
a[0][0] = a[1][1] = 1;
a[0][1] = a[1][0] = 0;
}
};
jz operator * (jz a, jz b) {
jz c;
fo(i, 0, 1) fo(j, 0, 1)
c.a[i][j] = (a.a[i][0] * b.a[0][j] + a.a[i][1] * b.a[1][j]) % mo;
return c;
}
jz operator + (jz a, jz b) {
fo(i, 0, 1) fo(j, 0, 1) a.a[i][j] = (a.a[i][j] + b.a[i][j]) % mo;
return a;
}
jz operator - (jz a, jz b) {
fo(i, 0, 1) fo(j, 0, 1) a.a[i][j] = (a.a[i][j] - b.a[i][j] + mo) % mo;
return a;
}
jz ksm(jz x, ll y) {
jz s;
for(; y; y /= 2, x = x * x)
if(y & 1) s = s * x;
return s;
}
jz z;
jz ksb(jz x, ll y) {
if(y == 1) return x;
if(y & 1) return ksb(x, y - 1) * x + x;
jz a = ksb(x, y / 2);
return a + a * ksm(x, y / 2);
}
jz calc(jz x, ll y) {
jz s = jz();
if(y > 0) s = s + ksb(x, y);
return s;
}
int main() {
scanf("%lld %d %d", &n, &k, &mo);
f[1] = f[2] = 1;
fo(i, 3, 3000000) {
f[i] = (f[i - 2] + f[i - 1]) % k;
if(!fi[f[i]]) fi[f[i]] = i;
}
work();
z.a[0][1] = z.a[1][0] = z.a[1][1] = 1; z.a[0][0] = 0;
ll ans = ksm(z, n).a[1][0];
if(l == -1) {
fo(i, 1, p0) if(p[i] <= n)
ans = (ans - ksm(z, n - p[i]).a[1][1] + mo) % mo;
} else {
fo(i, 1, l) if(p[i] <= n)
ans = (ans - ksm(z, n - p[i]).a[1][1] + mo) % mo;
if(p[l] <= n) {
ll v = (n - p[l]) / (p[r] - p[l]);
fo(i, l + 1, r) {
ll n2 = (ll) p[i] + v * (p[r] - p[l]);
if(n2 <= n) ans = (ans - ksm(z, n - n2).a[1][1] + mo) % mo;
}
if(v > 0) {
ll st = v * (p[r] - p[l]) + p[l];
jz sa; memset(sa.a, 0, sizeof sa.a);
fo(i, l + 1, r) sa = sa + ksm(z, p[r] - p[i]);
jz c = ksm(z, n - st) * sa * calc(ksm(z, p[r] - p[l]), v - 1);
ans = (ans - c.a[1][1] + mo) % mo;
}
}
}
pp("%lld
", ans);
}