Codeforces Round #566 (Div. 2) E. Product Oriented Recurrence

f_{x}=c^{2 x-6} \cdot f_{x-1} \cdot f_{x-2} \cdot f_{x-3} \text { for } x \geq 4

现在给出 $n, f_1, f_2, f_3, c,(4≤n≤10^{18}, 1≤f1, f2, f3, c≤10{9})$ 求 $f_n mod(10^9+ 7)$

先将式子变形

c^{x} f_{x}=c^{x-1} f_{x-1} \cdot c^{x-2} f_{x-2} \cdot c^{x-3} f_{x-3}

令 $g(x, p)=c^{x} f_{x}$ 中素数 $p$ 的个数，一个数可以看成若干个素数为基底的次幂形式，比如 $40=2^{3} \times 5$ （2的个数为3）。那么 $g(x, p)=g(x-1, p)+g(x-2, p)+g(x-3, p)$ （这实际上是一个取对数的形式），既然一个数可以由若干个素数组成，那么 $c^xf_x$ 可以表示成 $\sum _i {g(x,p_i)}$ 。

要得到 $c^nf_n$ 可以用矩阵快速幂完成，然后乘以 $c^n$ 的逆元就可以得到 $f_n$ 。

这里用矩阵快速幂有一个注意点，因为最终 $c^nf_n$ 结果表现为 $a_1^{c1}a_2^{c2}...$ ,这里矩阵快速幂用来计算素数基底的个数，也就是次幂项，就是 $c1,c2...$ ，对于其中任一项 $a^c$ 来说，结果要 $mod (10^9+7)$ ，由欧拉定理 $a^{\psi(p)}\equiv 1 (mod p)$ ，其中 $p$ 是素数；得 $a^{p-1}\equiv 1 (mod p)$ ，等式两边同时可以乘以若干个 $a^{p-1}$ ，即 $a^{k(p-1)}\equiv1(mod p)$ ，这时回到我们要求的 $a^c$ 上来，我们在前面这个等式两边同时乘以 $a^c$ ，即可以得到 $a^{k(p-1)+c} \equiv a^c (mod p)$ ，在题目中，我们通过快速幂可以得到的是 $k(p-1)+c$ ，将这个结果模 $p-1$ 就能得到在 $mod$ $p$ 下 $c$ 的大小。所以矩阵快速幂中要 $mod$ $(p-1)$ ，其余地方 $mod$ $p$ 。

另外在得到 $c^nf_n$ 后，通过费马小定理求 $c^n$ 的逆元时，我们 $mod$ $(p-2)$ 。

#include <bits/stdc++.h>
#define mem(a,b) memset(a, b, sizeof(a))
#define ll long long
#define pb push_back
#define fi first
#define se second
#define pi acos(-1)
#define eps 0.0001
using namespace std;
const int maxn = 1e6+ 5;
const int mod = 1e9+6;
const int mod1 = 1e9+7;
ll f1,f2,f3, c;
ll n;
//矩阵快速幂
struct matrix{
  ll val[16][16];
  matrix(){memset(val,0,sizeof(val));}
  matrix operator*(const matrix & rhs)const{
    matrix ret;
    for(int i = 0; i < 16; i++)
      for(int j = 0; j < 16; j++)
        for(int k = 0; k < 16; k++)
          ret.val[i][j] += val[i][k] * rhs.val[k][j]%mod,
          ret.val[i][j] %= mod;
    return ret;
  }
};
matrix pow_mod(const matrix& b, ll p){
  matrix ret, cur = b;
  for(int i = 0; i < 16; i++)
    ret.val[i][i] = 1;
  while(p){
    if(p & 1)
      ret = ret * cur;
    cur = cur * cur;
    p>>=1;
  }
  return ret;
}
//快速幂
ll pow_mod1(ll a, ll p){
    ll ans = 1;
    while(p){
        if(p & 1) ans = ans * a % mod1;
        a = a * a % mod1;
        p >>= 1;
    }
    return ans;
}
int main(){
    cin >> n >> f1 >> f2 >> f3 >> c;
    matrix p;
    //快速幂矩阵
    p.val[0][0] = p.val[1][0] = p.val[2][0] = 1;
    p.val[0][1] = p.val[1][2] = 1;
    p = pow_mod(p, n-3);
    map<ll, ll> mp1,mp2,mp3;//分别计算f_x-1,f_x-2,f_x-3，对应特定的素数的次幂项
    set<ll> Set;//收集所有素数基数
    


    ll t3 = c;
    for(ll i = 2; i <= sqrt(t3) && t3 > 1; i++){
        if(t3 % i == 0){
            Set.insert(i);
        }
        while(t3 % i == 0){
            t3 /= i;
            mp3[i]+=3;
            mp2[i]+=2;
            mp1[i]+=1;
        }
    }
    if(t3 > 1) {
        Set.insert(t3);
        mp3[t3] += 3;
        mp2[t3] += 2;
        mp1[t3] += 1;
    }


     t3 = f3;
    for(ll i = 2; i <= sqrt(t3) && t3 > 1; i++){
        if(t3 % i == 0){
            Set.insert(i);
     
        }
        while(t3 % i == 0){
            t3 /= i;
            mp3[i]++;
        }
    }
    if(t3 > 1) {mp3[t3] += 1;Set.insert(t3);}

   


    ll t2 = f2;
    for(ll i = 2; i <= sqrt(t2) && t2 > 1; i++){
        if(t2 % i == 0){
            Set.insert(i);
          
        }
        while(t2 % i == 0){
            t2 /= i;
            mp2[i]++;
        }
    }
    if(t2 > 1) {mp2[t2] += 1;Set.insert(t2);}




    ll t1 = f1;
    for(ll i = 2; i <= sqrt(t1) && t1 > 1; i++){
        if(t1 % i == 0){
            Set.insert(i);

        }
        while(t1 % i == 0){
            t1 /= i;
            mp1[i]++;
        }
    }
    if(t1 > 1) {mp1[t1] += 1;Set.insert(t1);}





    ll ans = 1LL;
    for(auto w : Set){
        matrix temp;
        temp.val[0][0] = mp3[w];
        temp.val[0][1] = mp2[w];
        temp.val[0][2] = mp1[w];
        temp = temp * p;

        ans = ans * pow_mod1(w,temp.val[0][0]) % mod1;
    }
    ll suf = pow_mod1(c, n);
    suf = pow_mod1(suf, mod1-2);
    ans = suf * ans % mod1;
    cout << ans << endl;

    return 0;
}