题解
- 拿到手, 想到是递推转矩阵快速幂, 但是我想的复杂了
刚开始想的是用a[n] * b[n]的递推构造矩阵+ 分块 + 快速幂, 不过这时候开的矩阵是10 * 10一组最大多O(size^3) * log n
size为矩阵阶数 - 最早思路:
aibi=ai−1bi−1AxBx+ai−1AxBy+bi−1AyBx+AyBy
所以:
⎡⎣⎢⎢⎢⎢⎢⎢AxBx0000AxByAx000AyBx0Bx00By101000101⎤⎦⎥⎥⎥⎥⎥⎥∗⎡⎣⎢⎢⎢⎢⎢⎢ai−1bi−1ai−1bi−1AyBy⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢aibiaibiAyBy⎤⎦⎥⎥⎥⎥⎥⎥
记B=⎡⎣⎢⎢⎢⎢⎢⎢AxBx0000AxByAx000AyBx0Bx00By101000101⎤⎦⎥⎥⎥⎥⎥⎥ Si=⎡⎣⎢⎢⎢⎢⎢⎢aibiaibiAyBy⎤⎦⎥⎥⎥⎥⎥⎥
B∗Si−1=Si,那么Sum(n)=∑n−1i=0Si,AOD(n)是Sum(n)的首项 -
Sum(n)=∑n−1i=0Si=∑n−1i=0Bi∗S0=(∑n−1i=0Bi)∗S0
令 Tn=∑n−1i=0Bi,T1=E
那么 Tn=E+B∗Tn−1
构造分块矩阵递推:
[TnE]=[Tn−1E]∗[BE0E]=[T1E]∗[BE0E]n−1
所以 Sum(n)=[EE]∗[BE0E]n−1∗S0
矩阵10*10, T了。。。 - 然后直接在递推时候考虑 AOD(n)
即: AOD(n) = AOD(n-1) + aibi,再次递推简化5 * 5, 细节处理A, 还是太弱了啊, 矩阵练了这么多。。
code
#include <iostream>
#include <cstring>
using namespace std;
typedef long long ll;
const ll mod = 1000000007;
ll n, a0, ax, ay, b0, bx, by;
struct Matrix{
ll mat[5][5];
void clr(){
memset(mat, 0, sizeof mat);
}
void E(){
clr();
for(int i = 0; i < 5; ++i) mat[i][i] = 1;
}
Matrix operator * (const Matrix &b){
Matrix res;
res.clr();
for(int i = 0; i < 5; ++i)
for(int j = 0; j < 5; ++j)
for(int k = 0; k < 5; ++k)
res.mat[i][j] = (res.mat[i][j] + mat[i][k] * b.mat[k][j]) % mod;
return res;
}
Matrix operator ^ (ll b){
Matrix res, tmp;
res.E();
memcpy(tmp.mat, mat, sizeof mat);
while(b){
if(b & 1) res = res * tmp;
tmp = tmp * tmp;
b >>= 1;
}
return res;
}
};
int main(){
while(cin >> n){
cin >> a0 >> ax >> ay >> b0 >> bx >> by;
if(n == 0) {
cout << "0" << endl;
continue;
}
a0 %= mod;
ax %= mod;
ay %= mod;
b0 %= mod;
bx %= mod;
by %= mod;
Matrix init, tmp;
init.clr();
tmp.clr();
init.mat[0][0] = a0 * b0 % mod;
init.mat[0][1] = a0 * b0 % mod;
init.mat[0][2] = a0;
init.mat[0][3] = b0;
init.mat[0][4] = 1;
/*tmp底数**/
tmp.mat[0][0] = 1;
tmp.mat[1][0] = ax * bx % mod;
tmp.mat[2][0] = ax * by % mod;
tmp.mat[3][0] = ay * bx % mod;
tmp.mat[4][0] = ay * by % mod;
tmp.mat[1][1] = ax * bx % mod;
tmp.mat[2][1] = ax * by % mod;
tmp.mat[3][1] = ay * bx % mod;
tmp.mat[4][1] = ay * by % mod;
tmp.mat[2][2] = ax;
tmp.mat[4][2] = ay;
tmp.mat[3][3] = bx;
tmp.mat[4][3] = by;
tmp.mat[4][4] = 1;
init = init * (tmp ^ (n - 1));
cout << init.mat[0][0] << endl;
}
return 0;
}
— 待修改