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

Sequence

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 310    Accepted Submission(s): 99

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.

 

Input

     The first line has a number,T,means testcase.

     Each testcase has 5 numbers,including n,a,b,c,p in a line.

    1≤T≤10,1≤n≤1018,1≤a,b,c≤109 , p is a prime number,and p≤109+7 .

 

Output

     Output one number for each case,which is fn mod p.

 

Sample Input

  
  

   
   1
5 3 3 3 233
  
  

 

Sample Output

  
  

   
   190
  
  

 

Source

BestCoder Round #80

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

题目大意：求公式的结果

题目分析：显然取对数，递推式一眼化成矩阵乘法，然后快速幂，因为矩阵快速幂的结果是a的幂，所以矩阵快速幂取模的时候根据费马小定理要模p-1，不过有可能a模p==0，这个时候答案显然是0，但是有可能指数取模的结果是0，这个时候a^0变成1就不对了

#include <cstdio>
#include <cstring>
#define ll long long
ll n, a, b, c, p;
struct matrix  
{  
    ll m[5][5];  
};  
  
matrix multiply(matrix x, matrix y, ll MOD)  
{  
    matrix ans;  
    memset(ans.m, 0, sizeof(ans.m));  
    for(int i = 1; i <= 3; i++)  
        for(int j = 1; j <= 3; j++)  
            if(x.m[i][j])  
                for(int k = 1; k <= 3; k++)  
                    ans.m[i][k] = (ans.m[i][k] + x.m[i][j] * y.m[j][k]) % MOD;  
    return ans;  
}  
  
matrix quickmod(matrix a, ll x, ll MOD)  
{  
    matrix ans;  
    memset(ans.m, 0, sizeof(ans.m));  
    for(int i = 1; i <= 3; i++)  
        ans.m[i][i] = 1;  
    while(x)  
    {  
        if(x & 1)  
            ans = multiply(ans, a, p - 1);  
        x >>= 1;  
        a = multiply(a, a, p - 1);  
    }  
    return ans;  
}  

ll qpow(ll x, ll nn)
{
    ll res = 1;
    while(nn)
    {
        if(nn & 1)
            res = (res * x) % p;
        x = (x * x) % p;
        nn >>= 1;
    }
    return res;
}

int main()
{
    int T;
    scanf("%d", &T);
    while(T --)
    {
        scanf("%I64d %I64d %I64d %I64d %I64d", &n, &a, &b, &c, &p);
        if(n == 1)
            printf("1\n");
        else if(n == 2)
            printf("%I64d\n", qpow(a, b));
        else
        {
            if(a % p == 0)
                printf("0\n");
            else
            {
                matrix ans;
                ans.m[1][1] = (b * c + b) % (p - 1);
                ans.m[1][2] = b;
                ans.m[1][3] = 1;
                matrix mx;
                mx.m[1][1] = c; mx.m[1][2] = 1; mx.m[1][3] = 0;
                mx.m[2][1] = 1; mx.m[2][2] = 0; mx.m[2][3] = 0;
                mx.m[3][1] = b; mx.m[3][2] = 0; mx.m[3][3] = 1;
                ans = multiply(ans, quickmod(mx, n - 3, p - 1), p - 1);
                printf("%I64d\n", qpow(a, ans.m[1][1]));
            }
        }
    }
}