Blackops

初心易得，始终难守

0%

2017年浙江工业大学大学生程序设计迎新赛热身赛 E 递推（矩阵快速幂）

发表于 更新于 分类于 ， ，

时间限制：C/C++ 2秒，其他语言4秒
空间限制：C/C++ 131072K，其他语言262144K
64bit IO Format: %lld

题目描述
JYM和XJ转眼就从小学上了高中。在学习递推的时候，JYM在纸上随手写了一个递推关系式：an=2an-1，a0=0。写完这个递推式，JYM拿给XJ看，XJ觉得太过简单，于是大笔一挥，在等式右边又加了一个式子，变成了这样：an=2an-1+n2。JYM看到这个式子，想要算几个项来看看，可是一算就发现这个数据量太大了，你能帮他解决这个问题吗？
输入描述:
输入数据有多组（不超过100组数据），每组数据包含一个整数N<=1018
输出描述:
一个整数X，表示递推式第n项的值。由于数字太大，因此结果对于1000000009取模后输出。
示例1
输入

0
1
2
3
输出

0
1
6
21

题目链接：E.递推
看范围懂做法系列……明显的矩阵快速幂啊，主要就是构造转移矩阵使得$(n-1)^2$变成$n^2$，通过YY+观察可以发现$n^2=(n-1)^2+2*(n-1)+1$，那么矩阵应该有三项，一个是平方数本身，一个是平方数下的底数，一个是$1$，然后就YY构造出了以下矩阵：

即

那么在$N>=2$的时候，我们可以用矩阵快速幂优化转移矩阵的乘法，中间过程取模用带取模的快速乘法即可。退役老年人果然只能做做新生赛了……
代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#include <bits/stdc++.h>
using namespace std;
#define INF 0x3f3f3f3f
#define LC(x) (x<<1)
#define RC(x) ((x<<1)+1)
#define MID(x,y) ((x+y)>>1)
#define fin(name) freopen(name,"r",stdin)
#define fout(name) freopen(name,"w",stdout)
#define CLR(arr,val) memset(arr,val,sizeof(arr))
#define FAST_IO ios::sync_with_stdio(false);cin.tie(0);
typedef pair<int, int> pii;
typedef long long LL;
const double PI = acos(-1.0);
const int N = 4;
const int mod = 1000000009;
inline int mul(int a, int b)
{
    int r = 0;
    while (b)
    {
        if (b & 1)
            r = (r + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return r;
}
struct Mat
{
    int A[N][N];
    void zero()
    {
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < N; ++j)
                A[i][j] = 0;
    }
    void one()
    {
        zero();
        for (int i = 0; i < N; ++i)
            A[i][i] = 1;
    }
    Mat operator*(const Mat &b)
    {
        Mat c;
        c.zero();
        for (int i = 0; i < N; ++i)
        {
            for (int k = 0; k < N; ++k)
            {
                if (A[i][k])
                {
                    for (int j = 0; j < N; ++j)
                    {
                        if (b.A[k][j])
                            c.A[i][j] = (c.A[i][j] + mul(A[i][k], b.A[k][j])) % mod;
                    }
                }
            }
        }
        return c;
    }
    friend Mat operator^(Mat a, LL b)
    {
        Mat r;
        r.one();
        while (b)
        {
            if (b & 1)
                r = r * a;
            a = a * a;
            b >>= 1;
        }
        return r;
    }
};
int main(void)
{
    LL n;
    int a[4] = {0, 1, 6, 21};
    while (~scanf("%lld", &n))
    {
        if (n <= 3)
            printf("%d\n", a[n]);
        else
        {
            Mat A, M;
            A.zero();
            M.zero();
            A.A[0][0] = a[1]; A.A[0][1] = 1; A.A[0][2] = 1; A.A[0][3] = 1;
            A.A[1][0] = a[0]; A.A[1][1] = 1; A.A[1][2] = 1; A.A[1][3] = 1;
            M.A[0][0] = 2;
            M.A[1][0] = 1; M.A[1][1] = 1;
            M.A[2][0] = 2; M.A[2][1] = 2; M.A[2][2] = 1;
            M.A[3][0] = 1; M.A[3][1] = 1; M.A[3][2] = 1; M.A[3][3] = 1;
            M = M ^ (n - 1);
            A = A * M;
            printf("%d\n", A.A[0][0]);
        }
    }
    return 0;
}