codeforces 1182E Product Oriented Recurrence 矩阵快速幂

题意：设f(n) = c ^ (2n - 6) * f(n - 1) * f(n - 2) * f(n - 3), 问第n项是多少？

思路：官方题解：我们先转化一下，令g(x) =  c ^ x * f(x), 那么原式转化为了g(x) = g(x - 1) * g(x - 2) * g(x - 3)。之后我们可以考虑把f(1), f(2), f(3)和c的质因子找出来，枚举质因子对答案的贡献。我们发现，如果是质因子的数目的话，乘法就变成了加法（相当于统计质因子的指数），这样就可以用矩阵乘法优化了。注意，矩阵转移的时候，模数是1e9 + 6，因为转移的时候是指数（欧拉定理）。其实基于这种想法，我们可以不用处理质因子，直接计算g(1), g(2)，g(3)对答案的贡献。

代码：

#include <bits/stdc++.h>
#define LL long long
using namespace std;
const LL mod = 1e9 + 7;
const LL mod1 = 1e9 + 6; 
map<LL, LL> mp[4];
set<LL> s;
set<LL>::iterator it; 
struct Matrix {
	LL a[3][3];
	
	void init(LL num = 0) {
		memset(a, 0, sizeof(a));
		for (int i = 0; i < 3; i++)
			a[i][i] = num;
	}
	
	Matrix operator * (const Matrix& now) const {
		Matrix res;
		res.init();
		for (int i = 0; i < 3; i++)
			for (int j = 0; j < 3; j++)
				for (int k = 0; k < 3; k++)
					res.a[i][j] = (res.a[i][j] + (a[i][k] * now.a[k][j]) % mod1) % mod1;
		return res; 
	}
	
	Matrix operator ^ (const LL num) const {
		Matrix ans, x = *this;
		ans.init(1);
		LL now = num;
		for (; now; now >>= 1) {
			if(now & 1) ans = ans * x;
			x = x * x;
		}
		return ans;
	}
	
	void print() {
		for (int i = 0; i < 3; i++) {
			for (int j = 0; j < 3; j++) {
				printf("%lld ", a[i][j]);
			}
			printf("\n");
		}
			
	} 
};
void div(LL num, int pos) {
	for (LL i = 2; i * i <= num; i++) {
		if(num % i == 0) {
			s.insert(i);
			while(num % i == 0) {
				mp[pos][i]++;
				num /= i;
			}
		}
	}
	if(num > 1) {
		s.insert(num);
		mp[pos][num]++;
	}
}
LL qpow(LL x, LL y) {
	LL ans = 1;
	for (; y; y >>= 1ll) {
		if(y & 1ll) ans = (ans * x) % mod;
		x = (x * x) % mod;
	}
	return ans;
}
LL a[4];
int main() {
	LL n;
	scanf("%lld", &n);
	for (int i = 0; i < 4; i++) {
		scanf("%lld", &a[i]);
		div(a[i], i);
	}
	Matrix x, y, x1;
	x.init();
	x.a[0][2] = x.a[1][2] = x.a[2][2] = x.a[2][1] = x.a[1][0] = 1;
	x = x ^ (n - 1);
	LL ans = 1;
	for (it = s.begin(); it != s.end(); it++) {
		y.init();
		for (int i = 0; i < 3; i++) {
			y.a[0][i] = mp[i][*it];
		}
		for (int i = 0; i < 3; i++)
			y.a[0][i] = (y.a[0][i] + ((LL)(i + 1) * mp[3][*it] % mod)) % mod;
		y = y * x;
		ans = (ans * qpow(*it, y.a[0][0]) % mod) % mod;
	}
	ans = ans * qpow(qpow(a[3], mod - 2), n) % mod;
	printf("%lld\n", ans);
}