Codeforces 653G Move by Prime 组合数学

题意：

有一个长度为\(n\)的正整数序列\(a\)，有这样一种操作：
每次可以选序列中的某一个数乘上或除以某一个素数。
求对于每一个子序列使其所有元素相等的最少操作次数之和。

分析：

因为两个素数之间互不影响，单独考虑每一个素数\(p\)。
设当前子序列的长度为\(k\)，对应的指数为\(e_1, e_2 \cdots e_k\)。
每次操作会将某一个\(e_i\)增加\(1\)或减少\(1\)。
将\(e_i\)对应到数轴上的点，每次操作就相当于让某个点向左或向右移动一个单位长度。
让它们都相等的最少操作次数就等于这些点到中位数的距离之和\(d\)。

• \(k=4\)时，\(d=e_4+e_3-e_2-e_1\)，中位数为\(e_2\)或\(e_3\)。
• \(k=5\)时，\(d=e_5+e_4-e_2-e_1\)，中位数为\(e_3\)。

下面考虑所有的子序列：
对于\(e_k\)，如果它在一个子序列的左半部分那么它对答案的贡献是\(-1\) ，如果它在序列的右半部分那么它对答案的贡献是\(+1\)。
可以在\(e_1, e_2 \cdots e_{k-1}\)中选若干数，以及在\(e_{k+1} \cdots e_n\)选若干数，构成包含\(e_k\)的子序列。

考虑下面这个生成函数：

\[(1+\frac{1}{x})^{k-1}+(1+x)^{n-k}=\frac{(1+x)^{n-1}}{x^{k-1}} \]

• \(e_k\)在左半部分的子序列的个数为指数为负的系数之和
• \(e_k\)在右半部分的子序列的个数为指数为正的系数之和

因此\(e_k\)对所有子序列的贡献为:

\[\sum\limits_{i=k}^{n-1}C_{n-1}^i-\sum\limits_{i=0}^{k-2}C_{n-1}^i \]

\[=\sum\limits_{i=0}^{n-1}C_{n-1}^i-\sum\limits_{i=0}^{k-1}C_{n-1}^i-\sum\limits_{i=0}^{k-2}C_{n-1}^i \]

\[=2^{n-1}-2\sum\limits_{i=0}^{k-2}C_{n-1}^i-C_{n-1}^{k-1} \]

\[=2^{n-1}-\sum\limits_{i=0}^{k-1}C_n^i \]

而且\({e_i}\)最大不会超过\(20\)，所以统计一下每个\(e_i\)出现的次数，再预处理一下组合数的前缀和的前缀和。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
const LL MOD = 1000000007;

LL mul(LL a, LL b) { return a * b % MOD; }

void add(LL& a, LL b) { a += b; if(a >= MOD) a -= MOD; }

void sub(LL& a, LL b) { a -= b; if(a < 0) a += MOD; }

LL pow_mod(LL a, int p) {
	LL ans = 1;
	while(p) {
		if(p & 1) ans = mul(ans, a);
		a = mul(a, a);
		p >>= 1;
	}
	return ans;
}

const int maxn = 300000 + 10;
const int maxp = 26000;

bool vis[maxn];
int prime[maxp], pid[maxn], pcnt;

void preprocess() {
	for(int i = 2; i < maxn; i++) {
		if(!vis[i]) { prime[pcnt] = i; pid[i] = pcnt++; }
		for(int j = 0; j < pcnt; j++) {
			if(i * prime[j] >= maxn) break;
			vis[i * prime[j]] = true;
			if(i % prime[j] == 0) break;
		}
	}
}

int n, a[maxn];
int cnt[maxp][20];

LL fac[maxn], inv[maxn], Cn[maxn];

void decompose(int x) {
	for(int i = 0; x > 1; i++) {
		int p = prime[i];
		if(p * p > x) break;
		if(x % p != 0) continue;
		int e = 0;
		while(x % p == 0) { x /= p; e++; }
		cnt[i][e]++;
	}
	if(x > 1) cnt[pid[x]][1]++;
}

int main()
{
	preprocess();

	scanf("%d", &n);
	for(int i = 1; i <= n; i++) scanf("%d", a + i);

	fac[0] = 1;
	for(int i = 1; i <= n; i++) fac[i] = mul(fac[i - 1], i);
	inv[n] = pow_mod(fac[n], MOD - 2);
	for(int i = n - 1; i >= 0; i--) inv[i] = mul(inv[i + 1], i + 1);
	for(int i = 0; i <= n; i++) Cn[i] = mul(mul(fac[n], inv[i]), inv[n - i]);
	for(int i = 1; i <= n; i++) add(Cn[i], Cn[i - 1]);
	for(int i = 1; i <= n; i++) add(Cn[i], Cn[i - 1]);

	for(int i = 1; i <= n; i++) decompose(a[i]);

	LL ans = 0;
	LL S = pow_mod(2, n - 1);
	for(int i = 0; i < pcnt; i++) {
		int tot = n;
		for(int j = 1; j < 20; j++) tot -= cnt[i][j];
		for(int j = 1; j < 20; j++) if(cnt[i][j]) {
			int L = tot, R = L + cnt[i][j] - 1;
			tot += cnt[i][j];
			LL t = Cn[R];
			if(L) sub(t, Cn[L - 1]);
			sub(t, mul(S, cnt[i][j]));
			add(ans, mul(t, j));
		}
	}

	printf("%lld\n", ans);

	return 0;
}