矩阵快速幂（模板）

转载请注明出处：http://blog.csdn.net/u012860063?viewmode=contents

Description

You have to find the nth term of the following function:

f(n) = a * f(n-1) + b * f(n-3) + c, if(n > 2)

= 0, if(n ≤ 2)

Input

Input starts with an integer T (≤ 100), denoting the number of test cases.

Each case contains four integers n (0 ≤ n ≤ 108),a b c (1 ≤ a, b, c ≤ 10000).

Output

For each case, print the case number and f(n) modulo 10007.

代码如下：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

struct Matrix
{
    int  m[5][5];
} I,A,T;
const int mod = 10007;
int a,b,n,mod;
Matrix matrixmul(Matrix a,Matrix b)
{
    int i,j,k;
    Matrix c;
    for (i=1; i<=4; i++)
        for(j=1; j<=4; j++)
        {
            c.m[i][j]=0;
            for(k=1; k<=4; k++)
            {
                c.m[i][j]+=(a.m[i][k]*b.m[k][j]);
                c.m[i][j]%=mod;
            }
        }
    return c;
}
Matrix quickpagow(int n)
{
    Matrix m=A,b=I;
    while (n>=1)
    {
        if(n&1)
            b=matrixmul(b,m);
        n=n>>1;
        m=matrixmul(m,m);
    }
    return b;
}

int main()
{
    int t,ans,i,flag=1,a,b,c;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d%d%d",&n,&a,&b,&c);
        memset(I.m,0,sizeof(I.m));
        memset(A.m,0,sizeof(A.m));
        for(i=1; i<=4; i++)
        {
            I.m[i][i]=1;
        }
        A.m[1][1]=a,A.m[1][2]=0,A.m[1][3]=b,A.m[1][4]=c;
        A.m[2][1]=A.m[3][2]=A.m[4][4]=1;
        if(n>=3)
        {
            T=quickpagow(n-2);//注意n-2为负的情况
            ans=0;
            printf("Case %d: %d\n",flag++,T.m[1][4]%mod);
        }
        else
            printf("Case %d: %d\n",flag++,0);
    }

    return 0;
}

Description

Let's define another number sequence, given by the following function:

f(0) = a

f(1) = b

f(n) = f(n-1) + f(n-2), n > 1

When a = 0 and b = 1, this sequence gives the Fibonacci sequence. Changing the values of a and b, you can get many different sequences. Given the values of a, b, you have to find the last m digits of f(n).

Input

Input starts with an integer T (≤ 10000), denoting the number of test cases.

Each test case consists of a single line containing four integers a b n m. The values of a and b range in [0,100], value of n ranges in [0, 109] and value of m ranges in [1, 4].

Output

For each case, print the case number and the last m digits of f(n). However, do NOT print any leading zero.

Sample Input

4

0 1 11 3

0 1 42 4

0 1 22 4

0 1 21 4

Sample Output

Case 1: 89

Case 2: 4296

Case 3: 7711

Case 4: 946

代码如下：

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
struct Matrix
{
	int mat[2][2];
};
Matrix A, f;
int n,mod;
Matrix mul(Matrix a,Matrix b)
{
	Matrix t;
	memset(t.mat,0,sizeof(t.mat));
	for(int i = 0; i < n; i++)
	{
		for(int k = 0; k < n; k++)
		{
			if(a.mat[i][k])
				for(int j = 0; j < n; j++)
				{
					t.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
					t.mat[i][j]%=mod;
				}
		}
	}
	return t;
}
Matrix quickP(int k)
{
	Matrix p = A ,m;
	memset(m.mat,0,sizeof(m.mat));
	for(int i=0;i<n;++i)//单位矩阵
	{
		m.mat[i][i]=1;
	}
	//if(k==0)
	//return e;
	while(k)
	{
		if(k & 1)
			m=mul(m,p);
		p=mul(p,p);
		k >>= 1 ;
	}
	return m;
}
int main()
{
	n = 2;//初值为2X2的矩阵
	int k,t;
	int a,b;//f[0]= a, f[1] = b;
	int z,w;
	int cas=0;
	while(scanf("%d",&t)!=EOF)
	{
		while(t--)
		{
			scanf("%d%d%d%d",&a,&b,&z,&w);
			mod=1;
			A.mat[1][1]=0;//初始化等比矩阵
			A.mat[0][0]=A.mat[1][0]=A.mat[0][1]=1;
			for(int i = 1 ; i <= w ; i++)//保留的位数
				mod*=10;
			Matrix ret=quickP(z-1);//z-1 乘的次数,注意z-1为负的情况
			int ans=(ret.mat[0][0]*b%mod+ret.mat[0][1]*a%mod)%mod;
			printf("Case %d: %d\n",++cas,ans);
		}
	}
	return 0;
}

----------------------------------------------------