[BZOJ3601]一个人的数论

[BZOJ3601]一个人的数论

试题描述

输入

见“试题描述”

输出

见“试题描述”

输入示例

见“试题描述”

输出示例

见“试题描述”

数据规模及约定

见“试题描述”

题解

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

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

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

\begin{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)\)

\begin{equation}
= \sum_{t|n} { \mu(t) t^d \sum_{i=1}^{d+1} { c_i \left( \frac{n}{t} \right) ^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\n", ans);
	
	return 0;
}