HDU-1005(矩阵二分快速幂)

来源：互联网发布：万网域名编辑：程序博客网时间：2024/05/22 13:06

Number Sequence

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 78704 Accepted Submission(s): 18378

Problem Description

A number sequence is defined as follows:

f(1) = 1, f(2) = 1, f(n) = (A * f(n - 1) + B * f(n - 2)) mod 7.

Given A, B, and n, you are to calculate the value of f(n).

Input

The input consists of multiple test cases. Each test case contains 3 integers A, B and n on a single line (1 <= A, B <= 1000, 1 <= n <= 100,000,000). Three zeros signal the end of input and this test case is not to be processed.

Output

For each test case, print the value of f(n) on a single line.

Sample Input

1 1 31 2 100 0 0

Sample Output

Author

CHEN, Shunbao

Source

ZJCPC2004

Recommend

JGShining

这道题目以前就已经做过,有印象,好像是49一循环,第一次交是按照7次一循环交上去的，错了，然后换了组数据打印出来一看，不行，要49一循环才行，但其实这都是你以前做过的原因，其实这样的题目是肯定有循环节的（都这么说）。。。

但是这是放在了快速幂里面的，怎么能是这样解题呢，自己是想不出来，所以看了一下别人的代码，才发现了用的矩阵，好神奇。。。其实每次自己得到一些新鲜的知识都是很神奇，神奇规神奇，但是你掌握不掌握就又是另一回事了。。。

自己讲的不如别人好，就不如把别人的讲解拿过来，给大家也和自己看一看

矩阵乘法与线性递推：

这里不再介绍矩阵的知识，相关可以参考线性代数。

还是从斐波纳契数列说起。前面说到，斐波纳契数列有一个通项公式：

F(n)=（√5/5)*{[(1+√5)/2]^n - [(1-√5)/2]^n}（√5表示根号5）。

对于上式的计算可以用前面说到的二分快速求幂发。现在我们推广到线性递推的快速幂计算。对于很多递推方程，很难直接计算出通项公式。这个时候如果计算很大的项的时候要借助矩阵乘法，然后对矩阵乘法进行二分求幂加速。而这种方法的关键是找出递推辅助矩阵：

设辅助矩阵为B 则应该满足下列条件

斐波纳契：{a(n-2) , a(n-1)} * B= {a(n-1) , a(n)} ;

容易得到B={0,1,1,1} (矩阵是2*2的，这里不方便打出分别为 X[1][1] X[1][2] X[2][1] X[2][2])

然后A = {1,1} 代表{a（1）， a(2)}

那么 A*B 就得到了 {a(2) , a(3)}

以此类推如果要得到第N项则是 result = A*B*B*..*B (n-2个B)。

然后根据矩阵乘法的结合律，先用二分求幂求出B*...*B部分然后再和A相乘在最后的结果矩阵result中result[1][2] = a(n);

这样就利用二分求幂和矩阵乘法得到了线性递推的快速解法。当然只要可以找到递推的辅助矩阵什么递归式都是可以转化成这种解法的。但是矩阵的维数太大的话也是会影响效率的。

经典题目：poj3070

别人的矩阵乘法模版：

#include<iostream>#include<cstdio>using namespace std;struct Matrix{int a[2][2];};Matrix E;void InitE(int size)       //初始化单位矩阵{for(int i=0;i<size;i++){for(int j=0;j<size;j++)E.a[i][j]=(i==j);}}Matrix MatrixMul(Matrix a,Matrix b)     //两矩阵相乘   {Matrix c;int i,j,k;for(i=0;i<2;i++){for(j=0;j<2;j++){c.a[i][j]=0;for(k=0;k<2;k++){c.a[i][j] += ((a.a[i][k]%10000)*(b.a[k][j]%10000));c.a[i][j]%=10000;}}}return c;}Matrix MatrixPow(Matrix a,int n)         //矩阵快速二分求n次幂{Matrix t=E;while(n>0){if(1&n)     //n是奇数t=MatrixMul(t,a);a=MatrixMul(a,a);n >>= 1;}return t;}int main(void){int a0,a1,p,q,k;Matrix matrix,m;InitE(2);while(scanf("%d %d %d %d %d",&a0,&a1,&p,&q,&k)!=EOF){if(!k)        //这两个特殊情况当时未考虑，导致多次WAprintf("%d\n",a0);else if(k==1)printf("%d\n",a1);else{matrix.a[0][0]=p;     //构造初始矩阵matrix.a[0][1]=q;matrix.a[1][0]=1;matrix.a[1][1]=0;m=MatrixPow(matrix,k-1);printf("%d\n",(m.a[0][0]*a1+m.a[0][1]*a0)%10000);}}return 0;}

自己的：

#include <stdio.h>#include <string.h>#include <iostream>#include <string>using namespace std;int digit[222];int A, B, N;struct Matrix{int m[2][2];}t;void build(int a, int b){t.m[0][0] = 0;t.m[0][1] = b % 7;t.m[1][0] = 1;t.m[1][1] = a % 7;}Matrix Mul(Matrix &a, Matrix &b){Matrix e;int ans = 0;for (int i = 0; i < 2; i++){for (int j = 0; j < 2; j++){for (int k = 0; k < 2; k++){ans += (a.m[i][k] % 7) * (b.m[k][j] % 7) % 7;}e.m[i][j] = ans % 7;ans = 0;}}return e;}Matrix binary_pow(Matrix temp, int n){int cnt = 0;Matrix e;while (n){digit[cnt++] = n % 2;n >>= 1;}e.m[0][0] = 1;e.m[0][1] = 0;e.m[1][0] = 0;e.m[1][1] = 1;for (int i = cnt - 1; i >= 0; i--){e = Mul(e, e);if (digit[i] == 1){e = Mul(e, temp);}}/*for (int i = 0; i < 2; i++){for (int j = 0; j < 2; j++){cout << e.m[i][j] << " ";}cout << endl;}*/return e;}int main(){while (scanf("%d%d%d", &A, &B, &N) != EOF){if (A == 0 && B == 0 && N == 0){break;}build(A, B);Matrix pow_ans = binary_pow(t, N - 1);int ans = (pow_ans.m[0][0] + pow_ans.m[1][0]) % 7;printf("%d\n", ans);}//system("pause");return 0;}