POJ3420 Quad Tiling

题意不多说，直接给出公式：

f(1)=1,f(2)=5,f(3)=11,f(4)=36

n≥5时 f(n)=f(n−1)+5∗f(n−2)+f(n−3)−f(n−4)

得到递推矩阵：

⎡⎣⎢⎢⎢⎢fn−3fn−2fn−1fn⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢000−1100101050014⎤⎦⎥⎥⎥∗⎡⎣⎢⎢⎢⎢fn−4fn−3fn−2fn−1⎤⎦⎥⎥⎥⎥

然后用矩阵快速幂的板子飞一飞就行了。。。

要注意的是，n≤4时，也需要对m取模，这个害我WA了半页。。。

#include <iostream>
#include <cstring>

using namespace std;

const int MAXN = 5;
const int MAXM = 5;
long long mod = 0;

struct Matrix {
int n, m;   //n列m行
long long a[MAXN][MAXM];

Matrix() { clear(); }

Matrix(int _n) {   //单位矩阵
clear();
n = m = _n;
for (int i = 0; i < _n; i++)
a[i][i] = 1;
}

Matrix(const Matrix &b) {
n = b.n;
m = b.m;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
a[i][j] = b.a[i][j];
}

void clear() {
n = m = 0;
memset(a, 0, sizeof(a));
}

Matrix operator+(const Matrix &b) const {
Matrix tmp;
tmp.n = n;
tmp.m = m;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
tmp.a[i][j] = a[i][j] + b.a[i][j];
}
return tmp;
}

Matrix operator-(const Matrix &b) const {
Matrix tmp;
tmp.n = n;
tmp.m = m;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
tmp.a[i][j] = a[i][j] - b.a[i][j];
return tmp;
}

Matrix operator*(const Matrix &b) const {
Matrix tmp;
tmp.clear();
tmp.n = b.n;
tmp.m = b.m;
for (int i = 0; i < n; ++i)
for (int j = 0; j < b.m; ++j) {
for (int k = 0; k < m; ++k)
tmp.a[i][j] = (tmp.a[i][j] + ((a[i][k] % mod) * (b.a[k][j] % mod)) % mod) % mod;
}
return tmp;
}
Matrix& operator%(const long long &mod)  {
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) {
a[i][j] = a[i][j]%mod;
}
return *this;
}

};

Matrix pow_mod(const Matrix &a, long long i, long long n) {
if (i == 0) return Matrix(4);
Matrix temp = pow_mod(a, i >> 1, n);
temp = temp * temp ;
if (i & 1) temp = temp * a ;
return temp;
}
//ans:(a^i)%n

long long n, m;
Matrix mCoefficient, mAnsCoefficient, mF1toF4, mFn;  //系数矩阵

void preWork() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
mF1toF4.clear();
mF1toF4.n = 4;
mF1toF4.m = 1;
mF1toF4.a[0][0] = 1;
mF1toF4.a[1][0] = 5;
mF1toF4.a[2][0] = 11;
mF1toF4.a[3][0] = 36;
mCoefficient.clear();
mCoefficient.n = 4;
mCoefficient.m = 4;
mCoefficient.a[0][1] = mCoefficient.a[1][2] = mCoefficient.a[2][3]
= mCoefficient.a[3][1] = mCoefficient.a[3][3] = 1;
mCoefficient.a[3][0] = -1;
mCoefficient.a[3][2] = 5;
}

int main() {
preWork();
while (cin >> n >> m , n) {
if (n > 4) {
mod = m;
mAnsCoefficient = pow_mod(mCoefficient, n-4, m);
mFn = mAnsCoefficient * mF1toF4;
mFn.a[3][0]=mFn.a[3][0]%m;
while (mFn.a[3][0]<0) mFn.a[3][0]=mFn.a[3][0]+m;
cout << mFn.a[3][0] << "\n";
} else {
cout << mF1toF4.a[n - 1][0]%m << "\n";
}
}
}

//f(n)=f(n-1)+5*f(n-2)+f(n-3)-f(n-4)