2020牛客寒假算法基础集训营1 J. 缪斯的影响力 （矩阵快速幂/费马小定理降幂）

https://ac.nowcoder.com/acm/problem/200658

f(n) = f(n-1) * f(n-2) * a^b ,f的第一项是x，第二项是y。

试着推出第三项是x·y·a^b，第四项是x·y2·a^2b，第五项是x²·y³·a^4b，第六项是x³y⁵a^7b

可以发现x的指数成1 0 1 1 2 3，y的指数0 1 1 2 3 5，a的指数是0 0 b 2b 4b 7b。

x和y的指数为斐波那契数列，a的指数规律为，除去系数b，其第n项前两项之和+1。

由于数据范围很大，所以可以用矩阵快速幂求出x y a的指数的第n项是多少。

x和y的乘法矩阵比较好构造

1 1

1 0

对于b的乘法矩阵，因为f（n）= f（n-1）+ f（n-2）+ 1

所以可以构造为：

1 1 1

1 0 0

0 0 1

因为最终是要求x·y·a的某次幂，且mod = 1e9+7是一个素数，所以这里再矩阵快速幂求解的过程中用费马小定理进行降幂操作

因为a^p-1≡1%p，所以a^b%p = a^b%(p-1)%p，在矩阵快速幂计算过程中矩阵元素相乘对mod-1取模而不是mod

求出指数后，再用快速幂求解即可

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9+7;
struct Matrix {
  ll a[3][3];
  Matrix() { memset(a, 0, sizeof a); }
  Matrix operator*(const Matrix &b) const {//重载矩阵乘法
    Matrix res;
    for (int i = 0; i <= 2; ++i)
      for (int j = 0; j <= 2; ++j)
        for (int k = 0; k <= 2; ++k)
          res.a[i][j] = (res.a[i][j] + (a[i][k] * b.a[k][j] )%(mod-1) ) % (mod-1);
    return res;
  }
}base1,base2,B,X,Y;
 
void init(){
    base1.a[0][0] = base1.a[0][1] = base1.a[1][0] = 1;
    base2.a[0][0] = base2.a[0][1] = base2.a[0][2] = base2.a[1][0] = base2.a[2][2] = 1;
    B.a[2][0] = X.a[1][0] = Y.a[0][0] = 1;
}
Matrix qpow(Matrix a,ll n) {
  Matrix res;
  res.a[0][0] = res.a[1][1] = res.a[2][2] = 1;
  while (n) {
    if (n & 1) res = res * a;
    a = a * a;
    n >>= 1;
  }
  return res;
}
ll power(ll a,ll n){
    ll res=1;
    a = a % mod;
    while(n){
        if(n&1)res=res*a%mod;
        n>>=1;
        a=a*a% mod;
    }
    return res;
}
 
int main(){
    ll n,x,y,a,b;
    cin>>n>>x>>y>>a>>b;
    if(n == 1) {
        cout<<x%mod;
        return 0;
    }
    if(n == 2){
        cout<<y%mod;
        return 0;
    }
    if(x%mod==0||y%mod==0||a%mod==0){cout<<0;return 0;}
    x%=mod,y%=mod;
    init();
    Matrix b1 = qpow(base1,n-2);
    Matrix b2 = qpow(base2,n-2);
    X = b1*X,Y = b1*Y,B = b2*B;
    ll Cx = X.a[0][0],Cy = Y.a[0][0],Cb = B.a[0][0]; 
    a = power(a%mod,b);
    ll ans = ((power(x,Cx)*power(y,Cy)%mod)*power(a,Cb)%mod);
    cout<<ans;
    return 0;
}