https://ac.nowcoder.com/acm/contest/3002/J

题意：f(1)=x,f(2)=y,f(i)=f(i−1)⋅f(i−2)⋅abf(1)=x,f(2)=y,f(i)=f(i-1)\cdot f(i-2)\cdot a^b，求 f(n),1≤n,x,y,a,b≤1012f(n),1\leq n,x,y,a,b\leq 10^{12}，对 10000000071000000007 取模。

矩阵快速幂，欧拉降幂

f(n)=xp⋅yq⋅(ab)mf(n)=x^p\cdot y^q\cdot (a^b)^m，其中 p,qp,q 分别为斐波那契数列的第 n−2,n−1n-2,n-1 项。mm 为递归式 f(i)=f(i−1)+f(i−2)+1,f(1)=1,f(2)=2f(i)=f(i-1)+f(i-2)+1,f(1)=1,f(2)=2 的第 n−2n-2 项。

由于项数特别大，因此需要像快速幂一样计算，复杂度为 O(log⁡nO(\log n)。

求斐波那契数列的递推式为

(f(i)f(i−1))=(1110)⋅(f(i−1)f(i−2))\begin{pmatrix} f(i) \\ f(i-1) \end{pmatrix} =\begin{pmatrix} 1&1 \\ 1&0 \\ \end{pmatrix} \cdot\begin{pmatrix} f(i-1) \\ f(i-2) \\ \end{pmatrix}

设 T=(1110)T=\begin{pmatrix}1&1 \\1&0\end{pmatrix}，Ai=(f(i)f(i−1))A_i=\begin{pmatrix}f(i)\\f(i-1)\end{pmatrix}，则 Ai=T⋅Ai−1A_i=T\cdot A_{i-1}，即 Ai=Ti−1⋅A1A_i=T^{i-1}\cdot A_1，这时对于 Ti−1T^{i-1} 用矩阵快速幂求解。

第二个递推式为 f(i)=f(i−1)+f(i−2)+Cf(i)=f(i-1)+f(i-2)+C，其中 C=1C=1。

(f(i)f(i−1)C)=(111100001)⋅(f(i−1)f(i−2)C)\begin{pmatrix} f(i) \\ f(i-1) \\ C \\ \end{pmatrix} =\begin{pmatrix} 1&1&1 \\ 1&0&0 \\ 0&0&1 \\ \end{pmatrix} \cdot\begin{pmatrix} f(i-1) \\ f(i-2) \\ C \\ \end{pmatrix}

求幂次的时候，还要同时取模，由费马定理，aϕ(p)%p=1a^{\phi (p)}\%p=1，所以要对 ϕ(p)\phi(p) 即 p−1p-1 取模。

调了半天一直卡90。因为如果x对mod取模后为0，而指数p刚好对mod-1取模后也为0，那么就会变成0的0次方，为1.但事实上需要的是0的mod-1次方，应该为0.这里快速幂是正确的，0的0次方确实应该为1，但是由于这是一个特例，所以一旦取模导致指数变化，变为0就会使得答案错误。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include <bits/stdc++.h>
using namespace std;
#define ll long long
typedef pair<int, int> pii;
const int maxn = 3e5 + 10;
const int INF = 0x3f3f3f3f;
const ll mod = 1e9 + 7;
ll n, x, y, a, b;
typedef vector<ll> vec;
typedef vector<vec> mat;
mat mul(const mat &a, const mat &b) {
    mat c(a.size(), vec(b[0].size()));
    for (int i = 0; i < (int)a.size(); i++) {
        for (int j = 0; j < (int)b[0].size(); j++) {
            for (int k = 0; k < (int)a[0].size(); k++) {
                c[i][j] += a[i][k] * b[k][j] % (mod - 1);
                c[i][j] %= (mod - 1);
            }
        }
    }
    return c;
}
mat pow(mat a, ll n) {
    mat res(a.size(), vec(a.size()));
    for (int i = 0; i < (int)a.size(); i++)
        res[i][i] = 1;
    while (n) {
        if (n & 1) res = mul(res, a);
        a = mul(a, a);
        n >>= 1;
    }
    return res;
}
mat res1(2, vec(1)), res2(3, vec(1));
void solve(ll n) {
    mat A(2, vec(2)), B(3, vec(3));
    A[0][0] = 1;
    A[0][1] = 1;
    A[1][0] = 1;
    A[1][1] = 0;
    res1[0][0] = res1[1][0] = 1;
    res1 = mul(pow(A, n - 1), res1);
 
    B[0][0] = B[0][1] = B[0][2] = 1;
    B[1][0] = 1; B[1][1] = B[1][2] = 0;
    B[2][0] = B[2][1] = 0; B[2][2] = 1;
    res2[0][0] = 2; res2[1][0] = 1; res2[2][0] = 1;
    res2 = mul(pow(B, n - 1), res2);
}
inline ll ksm(ll a, ll b) {
    ll ret = 1;
    a %= mod;
    while (b) {
        if (b & 1)ret = (ret*a) % mod;
        a = (a*a) % mod;
        b >>= 1;
    }
    return ret;
}
int main() {
    cin >> n >> x >> y >> a >> b;
    ll c = ksm(a, b);
    if (n == 1) { cout << x % mod << endl; return 0; }
    if (n == 2) { cout << y % mod << endl; return 0; }
    solve(n - 2);
    ll p = res1[1][0] + mod - 1, q = res1[0][0] + mod - 1, m = res2[1][0] + mod - 1;
    ll ans = ksm(x, p);
    ans = ans * ksm(y, q) % mod;
    ans = ans * ksm(c, m) % mod;
    cout << ans << endl;
    return 0;
}