hdu5895 Mathematician QSC（欧拉降幂+快速幂）

题目链接： hdu5895 ( Mathematician QSC )

此题先通过矩阵快速幂求解 \(g(n*y)\)，再根据欧拉降幂公式 \(x^{g(n*y)}\% (s+1)=x^{g(n*y)\%(s+1)+(s+1)}\% (s+1)\) 得出答案。

构造矩阵递推式的方式有多种。一种技巧性较强的方法是类似斐波那契数列用二维矩阵快速幂求出 \(f_{n}\) 和 \(f_{n+1}\) ，再根据公式可推出 \(g_n=f_nf_{n+1}/2\) （技巧性在于 \(g_n\) 的这个公式我推不出来）。

另一种方法是构造能直接求解 \(g_n\) 的矩阵递推式。

\[g_n=g_{n-1}+f_n^2 ···············(1) \\ f_{n+1}^2=4f_n^2+4f_{n-1}f_{n}+f_{n-1}^2···(2)\\ f_{n}^2=f_n^2······················(3)\\ f_nf_{n+1}=2f_{n}^2+f_{n-1}f_{n}·········(4) \]

首先写出 \((1)\) 式；发现右侧多了 \(f_n\) ，便在左侧补上 \(f_{n+1}\) ，写出 \((2)\) 式；右侧又多了 \(f_{n-1}f_n\) 和 \(f_{n-1}^2\) 便在左侧补上 \(f_{n}f_{n+1}\) 和 \(f_n^2\) ，写出 \((3)(4)\) 式；依次类推……直至左侧的所有状态可由其前一状态推出。

于是写出矩阵递推式

\[\begin{pmatrix}g_0\\ f_{1}^2\\ f_0^2\\ f_0f_{1} \end{pmatrix} =\begin{pmatrix}0\\ 1\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix}g_n\\ f_{n+1}^2\\ f_n^2\\ f_nf_{n+1} \end{pmatrix} =\begin{pmatrix}1&1&0&0\\ 0&4&1&4\\ 0&1&0&0\\ 0&2&0&1 \end{pmatrix} \begin{pmatrix}g_{n-1}\\ f_{n}^2\\ f_{n-1}^2\\ f_{n-1}f_{n} \end{pmatrix} (n>0) \]

\(AC\) 代码：

/**
 * hdu5895 Mathematician QSC
 * 矩阵快速幂，欧拉降幂公式
 */

#include <iostream>
#include <climits>
#include <cmath>
#include <iomanip>
#include <vector>
using namespace std;

typedef long long LL;

const int M = 4;
struct Matrix
{
    int row, col;
    LL m[M][M];
    Matrix(int row, int col)
    {
        this->row = row, this->col = col;
        for (int i = 0; i < row; ++i)
            for (int j = 0; j < col; ++j)
                m[i][j] = 0;
    }
    Matrix mul(const Matrix &t, const LL mod = LLONG_MAX)
    {
        int row = this->row, col = t.col, cnt = this->col;
        Matrix res(row, col);
        for (int i = 0; i < row; ++i)
            for (int j = 0; j < col; ++j)
                for (int k = 0; k < cnt; ++k)
                    res.m[i][j] += m[i][k]*t.m[k][j]%mod,
        			res.m[i][j] %= mod;
        return res;
    }
    Matrix pow(LL n, const LL mod = LLONG_MAX)
    {
        int sz = row;
        Matrix a = *this, t(sz, sz);
        for (int i = 0; i < sz; ++i) t.m[i][i] = 1;
        while (n) {
            if (n&1) t = t.mul(a, mod);
            a = a.mul(a, mod);
            n >>= 1;
        }
        return t;
    }
};

LL qpow(LL a, LL b, LL mod = LLONG_MAX)
{
    a %= mod;
    LL t = 1;
    while (b) {
        if (b&1) t = t*a%mod;
        a = a*a%mod;
        b >>= 1;
    }
    return t%mod;
}

LL phi(LL n)
{
    LL ans = n;
    for (LL i = 2; i*i <= n; ++i) {
        if (n%i) continue;
        ans = ans/i*(i-1);
        while (n%i == 0) n /= i;
    }
    if (n != 1) ans = ans / n*(n-1);
    return ans;
}

int main()
{
    Matrix A(4, 4), m0(4, 4);
    A.m[0][0] = A.m[0][1] = A.m[1][2] = A.m[2][1] = A.m[3][3] = 1;
    A.m[1][1] = A.m[1][3] = 4;
    A.m[3][1] = 2;
    m0.m[1][0] = 1;

    int T;
    cin >> T;
    LL n, y, x, s;
    while (T--) {
        cin >> n >> y >> x >> s;
        int p = phi(s+1);
        LL t;
        t = A.pow(n*y, p).mul(m0, p).m[0][0] + p;
        cout << qpow(x, t, s+1) << endl;
    }
    return 0;
}