[BZOJ3601]一个人的数论

[BZOJ3601]一个人的数论

试题描述

输入

见“试题描述”

输出

见“试题描述”

输入示例

见“试题描述”

输出示例

见“试题描述”

数据规模及约定

见“试题描述”

题解

如果题目给出的是某个数的质因数分解，思路就应该是构造积性函数。

它给出的函数无从知道是否是积性函数，并且即使它是积性函数我们也无法方便地求得 (f_d(p^a))，其中 (p) 是质数；这是因为互质的数把这个东西变得很难求。

那么不妨莫比乌斯反演一波

egin{equation}
sum_{i=1}^n { [gcd(i, n) = 1] cdot i^d } \
= sum_{i=1}^n { i^d sum_{t|i, t|n} mu(t) } \
= sum_{t|n} { mu(t) sum_{i=1}^{frac{n}{t}} (it)^d } \
= sum_{t|n} { mu(t) t^d sum_{i=1}^{frac{n}{t}} i^d }
end{equation}

(sum_{i=1}^n i^d) 这个叫做自然数幂求和，它能够表示成一个关于 (n) 的 (d+1) 次多项式（这个证明应该显然，我们把 (n^d + (n-1)^d + (n-2)^d + cdots + 2^d + 1^d) 展开就会得到一个关于 (n) 的多项式，显然 (n^d) 项在合并同类项之前有 (n) 项，所以它就是 (d+1) 次多项式）；现在，不妨令这个多项式的 (i) 次项系数为 (c_i)，即

[sum_{i=1}^{d+1} c_i n^i = sum_{i=1}^n i^d ]

这个 (c_i) 可以通过暴力求出前 (d+1) 项，然后高斯消元得出。（想拉格朗日插值的话也可以）

那么我们可以继续我们的反演了，接 ((1))

egin{equation}
= sum_{t|n} { mu(t) t^d sum_{i=1}^{d+1} { c_i left( frac{n}{t} ight) ^i } } \
= sum_{i=1}^{d+1} { c_i sum_{t|n} { mu(t) t^{d-i} n^i } } \
= sum_{i=1}^{d+1} { c_i n^i sum_{t|n} mu(t) t^{d-i} }
end{equation}

令 (H_t(n) = sum_{x|n} mu(x) x^t)，不难发现 (mu(x))、(x^t) 都是积性函数，所以 (mu(x) x^t) 是积性函数，那么 (H_t(n)) 就是 (mu(x) x^t) 与 (1) 狄利克雷卷积后的函数，也是积性的。

于是有

[H_t(n) = prod_{j=1}^w H_t(p_j^{alpha_j}) ]

然后我们展开一下 (H_t(p^a)) 看看是什么样子（(p) 是质数）

[H_t(p^a) = sum_{i=0}^a mu(p^i) p^{it} = 1 - p^t ]

这是因为 (mu(1) = 1, mu(p) = -1)。

接 ((2))，把刚推的这一坨带进去

[= sum_{i=1}^{d+1} { c_i n^i prod_{j=1}^w H_{d-i}(p_j^{alpha_j}) } ]

然后就可以暴力求了。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 1010
#define maxd 110
#define MOD 1000000007
#define LL long long

int d, w, pr[maxn], ex[maxn];

int Pow(int a, int b) {
	int ans = 1, t = a;
	while(b) {
		if(b & 1) ans = (LL)ans * t % MOD;
		t = (LL)t * t % MOD; b >>= 1;
	}
	return ans;
}

int A[maxd][maxd], coef[maxd];
void elim(int *a, int *b, int i, int n) {
	int rate = MOD - (LL)b[i] * Pow(a[i], MOD - 2) % MOD;
	rep(j, 1, n) {
		b[j] += (LL)rate * a[j] % MOD;
		if(b[j] >= MOD) b[j] -= MOD;
	}
	return ;
}
void gauss(int n) {
	rep(i, 1, n) {
		int j = i;
		for(; j <= n; j++) if(A[j][i]) break;
		swap(A[i], A[j]);
		rep(j, i + 1, n) if(A[j][i]) elim(A[i], A[j], i, n + 1);
	}
	dwn(i, n, 1) {
		int s = A[i][n+1];
		rep(j, i + 1, n) {
			s -= (LL)coef[j] * A[i][j] % MOD;
			if(s < 0) s += MOD;
		}
		coef[i] = (LL)s * Pow(A[i][i], MOD - 2) % MOD;
	}
	return ;
}

int H(int t, int p) {
	if(t < 0) return (1 + MOD - Pow(Pow(p, -t), MOD - 2)) % MOD;
	return (1 + MOD - Pow(p, t)) % MOD;
}

int main() {
	d = read(); w = read();
	rep(i, 1, w) pr[i] = read(), ex[i] = read(); // */
	
	int s = 0;
	rep(i, 1, d + 1) {
		s += Pow(i, d); if(s >= MOD) s -= MOD;
		A[i][d+2] = s;
		int prod = 1;
		rep(j, 1, d + 1) {
			prod = (LL)prod * i % MOD;
			A[i][j] = prod;
		}
	}
	gauss(d + 1);
	
	int pown = 1, ans = 0;
	rep(i, 1, d + 1) {
		int prod = 1;
		rep(j, 1, w) {
			pown = (LL)pown * Pow(pr[j], ex[j]) % MOD;
			prod = (LL)prod * H(d - i, pr[j]) % MOD;
		}
		ans += (LL)coef[i] * pown % MOD * prod % MOD;
		if(ans >= MOD) ans -= MOD;
	}
	
	printf("%d
", ans);
	
	return 0;
}