HDU-5667 Sequence (矩阵快速幂)

最新推荐文章于 2022-08-12 15:32:06 发布

原创最新推荐文章于 2022-08-12 15:32:06 发布 · 590 阅读

0 ·

CC 4.0 BY-SA版权

数学专栏收录该内容

38 篇文章

订阅专栏

本文介绍了一种使用矩阵快速幂方法解决Fibonacci数列问题的算法，详细解释了如何通过矩阵相乘和快速幂运算来求解特定位置的数列元素，并提供了实现代码。

题目：

http://acm.hdu.edu.cn/showproblem.php?pid=5667

题意：

fn=⎧⎩⎨⎪⎪1,ab,abfcn−1fn−2,n=1n=2otherwise

按这个表达式求出某一位的数字

思路：

公式整理一下如果只看指数的话就是f(n) = c*f(n-1)+f(n-2)+1,直接求出这个递推关系式就行了，矩阵快速幂直接求。

c 1 1 f(n-1)

1 0 0 * f(n-2)

0 0 1 1

如上两个矩阵相乘就是数列的下一位，用矩阵处理这类递推关系的方法很简单。

然后就得到了指数，注意因为这个也要取模，因为费马小定理a^(p-1)%p==1，所以只要将指数也就是矩阵相乘时对p-1取模，最后再将a^b对结果求快速幂，得到的就是答案了。

代码：

//kopyh
#include <bits/stdc++.h>
#define INF 0x3f3f3f3f
#define MOD 1000000007
#define N 11
using namespace std;
long long n,m,sum,res,flag;
struct Matrix
{
    long long ma[N][N];
    long long n,m;
};
Matrix matMul(Matrix m1, Matrix m2, long long mod)
{
    Matrix ans;
    memset(ans.ma,0,sizeof(ans.ma));
    for(int i=0;i<m1.n;i++)
        for(int j=0;j<m2.m;j++)
            for(int k=0;k<m1.m;k++)
                ans.ma[i][j]=(ans.ma[i][j]+m1.ma[i][k]*m2.ma[k][j]+mod)%mod;
    ans.n=m1.n; ans.m=m2.m;
    return ans;
}
Matrix matPow(Matrix m1, long long k, long long mod)
{
    Matrix ans;
    for(int i=0;i<m1.n;i++)
        for(int j=0;j<m1.m;j++)
            ans.ma[i][j] = (i==j);
    ans.n=ans.m=m1.n;
    while(k)
    {
        if(k&1)ans = matMul(ans,m1,mod);
        m1 = matMul(m1,m1,mod);
        k>>=1;
    }
    return ans;
}
long long power(long long x,long long k,long long mod)
{
	long long ans = 1;
	while(k)
    {
		if(k & 1) ans=ans*x%mod;
		x=x*x%mod;
		k >>= 1;
	}
	return ans;
}
int main()
{
    long long i,j,k,cas,T,t,z,a,b,c,p;
    scanf("%I64d",&T);
    cas=0;
    while(T--)
    {
        scanf("%I64d%I64d%I64d%I64d%I64d",&n,&a,&b,&c,&p);
        Matrix x,y;
        x.ma[0][0]=c;
        x.ma[0][1]=1;
        x.ma[0][2]=1;
        x.ma[1][0]=1;
        x.ma[1][1]=0;
        x.ma[1][2]=0;
        x.ma[2][0]=0;
        x.ma[2][1]=0;
        x.ma[2][2]=1;
        x.n=x.m=3;
        y.ma[0][0]=1;
        y.ma[1][0]=0;
        y.ma[2][0]=1;
        y.n=3;y.m=1;
        x = matPow(x,n-2,p-1);
        x = matMul(x,y,p-1);
        sum = power(a,b,p);
        res = power(sum,x.ma[0][0],p);
        printf("%I64d\n",res%p);
    }
    return 0;
}