hdu 5667 Sequence(矩阵快速幂+费马小定理+快速幂)

Problem Description

     Holion August will eat every thing he has found.

     Now there are many foods,but he does not want to eat all of them at once,so he find a sequence.

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

     He gives you 5 numbers n,a,b,c,p,and he will eat  fn  foods.But there are only p foods,so you should tell him  fn  mod p.

题意：

 

问f(n)对p取模的结果;

 

思路：

 

( a ^b) ^p[n]= a ^b * ((a ^b) ^p[n-1]) ^c * ((a ^b) ^p[n-2]);递推式子可以这样写;

合并后变为(a ^b) ^p[n]=(a ^b) ^{(c*p[n-1]+p[n-2]+1);}

可以得到p[n]=c*p[n-1]+p[n-2]+1;

这样的递推式可以用矩阵乘法得到第n项;这是矩阵乘法的一个应用,给matrix67大神的博客地址可以学习,

构造矩阵乘法:

 

p[n]　　　　c    1    1　　　　p[n-1]

p[n-1]    =  1   0    0   *  　p[n-2]

1　　　　　  0   0    1　　　　　　1

 

然后这中间还有一个问题,就是取模的问题; a^b*p[n]%mod=a^{b*p[n]%(mod-1)}%mod;

这是根据费马小定理得到的;a ^(p-1)Ξ1%p;

 

a^b*p[n]%mod=a^{b*p[n]/(mod-1)*(mod-1)+b*p[n]%(mod-1)}%mod;

 

令x=b*p[n]/(mod-1)则a^b*p[n]%mod=a^{x*(mod-1)+b*p[n]%(mod-1)}%mod=a^{b*p[n]%(mod-1)}%mod;

这就是取模的结果啦;

然后就是快速幂算答案啦;

还有要注意的就是a%p==0的情况;

#include <cstdio>
#include <iostream>
#include <cstring>
#include <string>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <vector>
#include <set>
#include <list>
#include <queue>
#include <map>
#include <stack>
using namespace std;
#define L(i) i<<1
#define R(i) i<<1|1
#define INF  0x3f3f3f3f
#define pi acos(-1.0)
#define eps 1e-9
#define maxn 100010
#define MOD 1000000007
int N;
long long p,a,b,c;
struct Matrix
{
    long long a[22][22];
    Matrix()
    {
        memset(a,0,sizeof(a));
        for(int i = 0; i < N; i++)
            a[i][i] = 1;
    }
};
Matrix operator*(Matrix A,Matrix B)
{
    Matrix ans;
    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++)
    {
        ans.a[i][j] = 0;
        for(int k = 0; k < N; k++)
            ans.a[i][j] += (A.a[i][k] * B.a[k][j])%(p-1);
        ans.a[i][j] %= (p-1);
    }
    return ans;
}
Matrix Pow(Matrix A,long long n)
{
    Matrix ans;
    while(n)
    {
        if(n % 2)
            ans = ans * A;
        n /= 2;
        A = A * A;
    }
    return ans;
}
long long mod_pow(long long aa,long long bb,long long mod)
{
    long long ans = 1;
    while(bb)
    {
        if(bb&1)
            ans *= aa;
        ans %= mod;
        bb /= 2;
        aa *= aa;
        aa %= mod;
    }
    return ans;
}
int main()
{
	int t;
	scanf("%d",&t);
	while(t--)
	{
	    N = 3;
		long long n;
		scanf("%lld%lld%lld%lld%lld",&n,&a,&b,&c,&p);
		if(n==1)
        {
            printf("1\n");
            continue;
        }
        if(n==2)
        {
            printf("%lld\n",mod_pow(a,b,p));
            continue;
        }
		Matrix A;
		memset(A.a,0,sizeof(A.a));
		A.a[0][0] = c;
		A.a[0][1] = 1;
		A.a[0][2] = b;
		A.a[1][0] = 1;
		A.a[2][2] = 1;
		A = Pow(A,n-2);
		long long k = A.a[0][0] * b + A.a[0][2];
		k = mod_pow(a,k,p);
		printf("%lld\n",k);
	}
}