hdu5895 Mathematician QSC（矩阵快速幂+逆元+降幂）

hdu5895

题目

求题目里的一个式子的值

思路

首先f(n)是可以通过矩阵递推得到的
其次，g(n)也可以归纳得出g(n)=f(n)*f(n+1)/2
对于x^(g(n*y))%(s+1),幂的次数太高，我们要用降幂公式
x^k%p=x^(k%phi(p)+phi(p))%p,其中phi是欧拉函数
但是对于这题，g(n)的表达式中有/2,所以我们需要求逆元，
但这个式子可能没有逆元，那么怎么做呢，我们就用到了一个逆元的应用，
(a/b)%c=(a%(bc))/b,把这些公式综合一下，结合一下代码，然后就可以求解了。

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <queue>
#include <algorithm>
#include <set>

using namespace std;

typedef long long ll;

ll n,y,x,s;

struct Matrix
{
    ll mat[2][2];
};

ll Phi(ll x)
{
    ll ans = x;
    for(ll i=2ll; i*i<=x; i++)
    {
        if(x % i == 0)
        {
            ans -= ans/i;
            while(x % i == 0)
                x /= i;
        }
    }
    if(x > 1)
        ans -= ans/x;
    return ans;
}

Matrix P = {0, 1,
            1, 2,
           };

Matrix I = {1, 0,
            0, 1,
           };

Matrix Matrix_Mul(Matrix a, Matrix b, ll MOD)
{
    Matrix c;
    for(int i=0; i<2; i++)
        for(int j=0; j<i; j++)
            swap(b.mat[i][j], b.mat[j][i]);
    for(int i=0; i<2; i++)
    {
        for(int j=0; j<2; j++)
        {
            c.mat[i][j] = 0;
            for(int k=0; k<2; k++)
            {
                c.mat[i][j] += a.mat[i][k] * b.mat[j][k];
                c.mat[i][j] %= MOD;
            }
        }
    }
    return c;
}

Matrix Matrix_Quick_Mod(Matrix a, ll m, ll MOD)
{
    Matrix ans = I;
    while(m)
    {
        if(m & 1)
            ans = Matrix_Mul(ans, a, MOD);
        m>>=1;
        a = Matrix_Mul(a, a, MOD);
    }
    return ans;
}

ll Quick_Mod(ll a, ll b, ll MOD)
{
    ll ans = 1ll;
    while(b)
    {
        if(b & 1)
            ans = (ans*a) % MOD;
        b>>=1;
        a = (a*a) % MOD;
    }
    return ans;
}

int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%I64d%I64d%I64d%I64d",&n,&y,&x,&s);
        s+=1;
        ll mod=Phi(s);
        mod*=2ll;
        Matrix tmp1 = Matrix_Quick_Mod(P, n*y, mod);
        ll ans1 = ( ( (tmp1.mat[1][1]*tmp1.mat[1][0]) % mod ) >> 1ll );
        ans1+=mod;
        ll ans=Quick_Mod(x,ans1,s);
        printf("%I64d\n",ans);

    }
    return 0;
}