hdu2971 递推+矩阵快速幂

解题报告

题目 ：http://acm.hdu.edu.cn/showproblem.php?pid=2971

题目描述 ：a[n] = 2 * p * a[n – 1] – a[n – 2], 求 s[n] = ∑a[i]^2. (n = 10^9)

算法 ：递推+矩阵快速幂

看到是递推的式子，n又是10^9,肯定是构造关于an或者是sn的矩阵向量，然后构造辅助矩阵，利用矩阵快速幂解决。

下面来递推：

è-S[n] = a[1] *a[1] + a[2]*a[2] + a[3]*a[3] + …

è-S[n] = a[1]^2 + a[2]^2 + (2 * p * a[2] – a[1]) * (2 * p * a[2] – a[1]) * …….

è-S[n] = a[1]^2 + a[2]^2 + 4p^2∑a[n-1]^2 – 4p(a[n-1]*a[n-2]) + ∑a[n-2]^2

è-S[n] = a[1]^2 + a[2]^2 + 4p^2(S[n-1] – a[1]^2) – 4p(a[n-1]*a[n-2]) + S[n – 2]

è>令 t[n] = a[n] * a[n – 1], C = a[1]^2 + a[2]^2 – 4p^2*a[1]^2 (常数), T[n] = ∑t[n]

è-S[n] = C + 4p^2*S[n-1] – 4p*T[n] + S[n – 2]

è-这就得到了S[n]与S[n-1]，T[n]，S[n-2]之间的关系，下面找T[n]表达式

è>t[n] = a[n] * a[n – 1] = (2 * a[n-1] – a[n-2]) * a[n-1]

è>T[n] = ∑2 * a[n-1]^2 – t[n-1] (n >=2)

è>T[n] = 2 *(S[n-1] – a[1]^2) – T[n-1] + t[1](因为上边是n>=2，所以要加个t[1] = a[1])

è>T[n] = 2*S[n-1] – T[n-1] +t[1] – 2*a[1] =  2*S[n-1] – T[n-1] - a[1]

è>这样我们得到了T[n] 与 S[n-1],T[n-2]和某个常数之间的关系

è>综合S[n]，T[n]

è>【S[n-1],s[n-2],T[n-1],1】 * (辅助矩阵tem) = 【S[n], S[n-1,T[n], 1】

è.利用S[n] = C + 4p^2*S[n-1] – 4p*T[n] + S[n – 2]和T[n] = 2*S[n-1] – T[n-1] - a[1]填辅助矩阵

è>利用矩阵快速幂求tem的n-2次方，然后求解。

值得注意的地方是由于数据规模比较大，而mod运算比较慢，因此尽量少用mod，利用I64的大数据范围特点，我们将矩阵的行和列全部加完后在取模，能减省很多时间，如

for(i = 0; i < p; i ++)

       for(j = 0; j < r; j ++){

           for(k = 0; k < q; k ++)

              ans.row[i].line[j] += a.row[i].line[k] * b.row[k].line[j];

           ans.row[i].line[j] %= m;

       }

而不是:

for(i = 0; i < p; i ++)

       for(j = 0; j < r; j ++){

           for(k = 0; k < q; k ++)

              ans.row[i].line[j] =（ ans.row[i].line[j]+ a.row[i].line[k] * b.row[k].line[j]）%= m;

       }

提交情况 ：wrong answer 2次（辅助矩阵构造错误）

                      Time limit error 5次（上面Mod用的太多）

                      Accepted 1次

经验与收获

1：知道了mod运算的速度

2：矩阵乘法是可以优化的

3：递推式中具有常数是怎么处理

4：两个递推之间的结合

 

AC code：

#include <stdio.h>

#define MAXN 5

typedef __int64 I64;

 

struct LINE{

    I64 line[MAXN];

};

 

struct ROW{

    LINE row[MAXN];

};

ROW A, tem, B;

I64 m;

ROW Matrix_multi(ROW a, ROW b, int p, int q, int r){

    int i, j, k;

    ROW ans = {0};

    for(i = 0; i < p; i ++)

       for(j = 0; j < r; j ++){

           for(k = 0; k < q; k ++)

              ans.row[i].line[j] += a.row[i].line[k] * b.row[k].line[j];

           ans.row[i].line[j] %= m;

       }

    return ans;

}

 

 

ROW MOD(ROW temp, int n){

    ROW ans = {0};

    for(int i = 0; i < 4; i ++) ans.row[i].line[i] = 1;

    while(n){

       if(n & 1) ans = Matrix_multi(ans, temp, 4, 4, 4);

       temp = Matrix_multi(temp, temp, 4, 4, 4);

       n >>= 1;

    }

    return ans;

}

 

int main(){

    int n, T;

    I64 p;

    scanf("%d", &T);

    while(T --){

       scanf("%I64d %d %I64d", &p, &n, &m);

       A.row[0].line[0] = (1 + p * p) % m;

       A.row[0].line[1] = 1;

       A.row[0].line[2] = p % m;

       A.row[0].line[3] = 1;

       tem.row[0].line[0] = 4 * p * p % m; tem.row[0].line[2] = 2 * p % m;

       tem.row[1].line[0] = tem.row[0].line[1] = tem.row[3].line[3] = 1;

       tem.row[0].line[3] = tem.row[1].line[1] = tem.row[1].line[2] = tem.row[1].line[3] = 0;

       tem.row[2].line[1] = tem.row[2].line[3] = tem.row[3].line[1] = 0;

       tem.row[2].line[0] = (-4 * p % m + m) % m; tem.row[2].line[2] = m - 1;

       tem.row[3].line[0] = ((1 - 3 * p * p) % m + m) % m; tem.row[3].line[2] = (-p % m + m) % m;

       tem = MOD(tem, n - 2);

       if(n == 2) B = A;

       else B = Matrix_multi(A, tem, 1, 4, 4);

       printf("%I64d\n", B.row[0].line[0] % m);

    }

    return 0;

}