XDU-1111 Too Easy（矩阵快速幂）

1111: Too Easy

Time Limit: 1 Sec   Memory Limit: 128 MB
[ Submit][ Status][ Web Board]

Description

有n个白色正方形排成一排，现在你需要把至少2个正方形涂黑，并且相邻两个正方形不能都是黑色，求方案数（不考虑旋转和翻转），并对1000000007（1e9+7）取模。

Input

多组测试数据（大约10000组），处理到EOF。
每组数据包含一行，表示n（2  <= n <= 109）。

Output

对于每组数据输出一行，表示方案数对1000000007（1e9+7）取模。

Sample Input

2
3
4

Sample Output

0
1
3

做了一天，我真是太渣了。

先算了前几个，在oeis查到了这个数列（http://oeis.org/A001924），看到了下方的公式有：a[n]=Fib(n+4)-n-3，想到斐波那契数列有通项公式，就直接用double求出来，看到前几十个都是符合的就交了，结果WA，重新测试很大的数（44以后），由于有n次方运算，double的精度产生的误差就会放大到比1000000007大，从而影响运算结果。

在写完之前，经人点播，其实已经发现了递推公式，但是看到有递推公式中有n和常数，就放弃了矩阵快速幂（果真是做题太少了），更坚定的通过通项求解的信念。。。

其实在我自己数的时候就已经用了求递推公式的方法，只不是没注意到。。。

设g(x)=1+2+…+x；

刚开始只看到每次只求一个的个数，没有注意到他的个数与前面的个数的关系

f(3)=g(1)

f(4)=g(2)

f(5)=g(3)+g(1)=g(2)+3+g(1)=f(4)+f(3)+5

f(6)=g(4)+g(1)+g(2)=g(3)+4+g(1)+g(2)=g(3)+g(1)+g(2)+4=f(5)+f(4)+4

f(7)=g(5)+g(3)+g(1)+g(2)+g(1)=g(4)+5+g(3)+g(1)+g(2)+g(1)=g(4)+g(1)+g(2)+g(3)+g(1)+5=f(6)+f(5)+5

…………

可得递推公式为：f(n)=f(n-1)+f(n-2)+n-2;        ①

                                f(n-1)=f(n-2)+f(n-3)+n-3;     ②

①-②得：f(n)=2*f(n-1)-f(n-3)+1;               ③

                f(n-1)=2*f(n-2)-f(n-4)+1;             ④

③-④得：f(n)=3*f(n-1)-2*f(n-2)-f(n-3)+f(n-4)

则可以构造：

#include <cstdio>

using namespace std;

const int MAX=4;
const long long MOD=1000000007;

struct Matrix{
    long long  m[MAX][MAX];
};

Matrix P = {3,1,0,0,//此处由于编译器负数取模还是原数，所以先人工取模变成正数
            MOD-2,0,1,0,
            MOD-1,0,0,1,
            1,0,0,0};

Matrix I = {1,0,0,0,
            0,1,0,0,
            0,0,1,0,
            0,0,0,1};

Matrix matrixmul(Matrix a,Matrix b) {//矩阵乘法
    int i,j,k;
    Matrix c;
    for (i = 0 ; i < MAX; i++)
        for (j = 0; j < MAX;j++) {
            c.m[i][j] = 0;
            for (k = 0; k < MAX; k++)
                c.m[i][j] += (a.m[i][k] * b.m[k][j])%MOD;
            c.m[i][j] %= MOD;
        }
    return c;
}

Matrix quickpow(long long n) {
    Matrix m = P, b = I;
    while (n >= 1) {
        if (n & 1)
            b = matrixmul(b,m);
        n = n >> 1;
        m = matrixmul(m,m);
    }
    return b;
}

int main() {
    int n;
    long long f[7]={0,0,0,1,3,7,14};//f[n]分别表示n取3,4,5,6时的方案数
    Matrix ans;
    while(1==scanf("%d",&n)) {
        if(n<7)
            printf("%lld\n",f[n]);
        else {
            ans=quickpow(n-6);//用矩阵快速幂算递推结果
            printf("%lld\n",(f[6]*ans.m[0][0]+f[5]*ans.m[1][0]+f[4]*ans.m[2][0]+f[3]*ans.m[3][0])%MOD);
        }
    }
    return 0;
}