应该算一道结论题吧。
首先很明显的一点，就是对于ai,j，我们只用求出ai,1就行了，剩下的就是斐波那契和矩阵快速幂的事。
至于如何求ai,1，用两种方法，但是哪一个我都不懂为啥，网上也没有解释。
F1：斐波那契最小拆分。对于 i，拆分 i - 1：每一次选择小于i - 1且最大的斐波那契数相减，然后用二进制表示每一个斐波那契数取没取到。于是的到一个01序列，然后再后面再加一个01，对应到新的取数序列，代入斐波那契，就得到了ai,1。
举个例子：i = 5，于是拆分4：4 = 3 + 1，则序列为101，加上01：10101，对应到斐波那契序列为1 + 3 + 8 = 12，就是ai,1。
F2：解法就更迷了：把ai,1序列输到OEIS上，直接得到了这么个公式：ai,1 = i * (1 + √5) / 2 + i - 1。代入求值……
知道了ai,1和ai,2，剩下的就是矩阵快速幂水过了。
1 #include<cstdio> 2 #include<iostream> 3 #include<cmath> 4 #include<algorithm> 5 #include<cstring> 6 #include<cstdlib> 7 #include<cctype> 8 #include<vector> 9 #include<stack> 10 #include<queue> 11 using namespace std; 12 #define enter puts("") 13 #define space putchar(' ') 14 #define Mem(a, x) memset(a, x, sizeof(a)) 15 #define rg register 16 typedef long long ll; 17 typedef double db; 18 const int INF = 0x3f3f3f3f; 19 const db eps = 1e-8; 20 //const int maxn = ; 21 inline ll read() 22 { 23 ll ans = 0; 24 char ch = getchar(), last = ' '; 25 while(!isdigit(ch)) {last = ch; ch = getchar();} 26 while(isdigit(ch)) {ans = ans * 10 + ch - '0'; ch = getchar();} 27 if(last == '-') ans = -ans; 28 return ans; 29 } 30 inline void write(ll x) 31 { 32 if(x < 0) x = -x, putchar('-'); 33 if(x >= 10) write(x / 10); 34 putchar(x % 10 + '0'); 35 } 36 37 ll n, m, mod; 38 ll c[3]; 39 const int ha[] = {0, 3, 2, 3, 3, 2}; 40 41 const int N = 2; 42 struct Mat 43 { 44 ll a[N][N]; 45 Mat operator * (const Mat& oth)const 46 { 47 Mat ret; Mem(ret.a, 0); 48 for(int i = 0; i < N; ++i) 49 for(int j = 0; j < N; ++j) 50 for(int k = 0; k < N; ++k) 51 ret.a[i][j] += a[i][k] * oth.a[k][j], ret.a[i][j] %= mod; 52 return ret; 53 } 54 }f; 55 56 void init() 57 { 58 Mem(f.a, 0); 59 f.a[0][0] = f.a[0][1] = f.a[1][0] = 1; 60 } 61 62 Mat quickpow(Mat A, ll b) 63 { 64 Mat ret; Mem(ret.a, 0); 65 ret.a[0][0] = ret.a[1][1] = 1; 66 for(; b; b >>= 1, A = A * A) 67 if(b & 1) ret = ret * A; 68 return ret; 69 } 70 71 int main() 72 { 73 n = read(); m = read(); mod = read(); 74 c[1] = ((ll)(n * (1 + sqrt(5))) / 2 + n - 1) % mod; 75 c[2] = ((c[1] * 2 - n + 1) % mod + mod) % mod; 76 init(); 77 if(m < 3) {write(c[m]), enter; return 0;} 78 m -= 2; 79 Mat A = quickpow(f, m); 80 ll ans = A.a[0][0] * c[2] % mod + A.a[0][1] * c[1] % mod; 81 write(ans % mod), enter; 82 return 0; 83 }