『题解』Luogu-P6271 [湖北省队互测2014]一个人的数论

P6271 [湖北省队互测2014]一个人的数论

Description

• 给出非负整数 $d$ 和正整数 $n$ 的质因数分解式 $n={\prod }_{i=1}^{\omega }{p}_i^{{\alpha }_i}$，请求出
  $$\left(f_d(n)=\sum _{i=1}^n[\gcd (i,n)=1]i^d\right)\mathrm{m}\mathrm{o}\mathrm{d}(10^9+7)$$

• $0\le d\le 100,1\le \omega \le 10^3,2\le p_i\le 10^9,1\le {\alpha }_i\le 10^9$。

Solution

$$\begin{array}{rl}{f}_d(n)& =\sum _{i=1}^n[\gcd (i,n)=1]i^d\\ & =\sum _{i=1}^n{i}^d\sum _{k\mid \gcd (i,n)}\mu (k)\\ & =\sum _{k\mid n}\mu (k)\sum _{i=1}^n[k\mid i]i^d\\ & =\sum _{k\mid n}\mu (k)k^d\sum _{i=1}^{\frac{n}{k}}{i}^d\end{array}$$

后面的等幂和直接用伯努利数做就可以了。

我们知道
$$\begin{array}{rl}{S}_m(n)& =\sum _{k=0}^{n-1}{k}^m\\ & =\frac{1}{m+1}\sum _{k=0}^m\left(\genfrac{}{}{0ex}{}{m+1}{k}\right)B_k{n}^{m+1-k}\end{array}$$
设出系数
$$f_i=\frac{1}{d+1}\left(\genfrac{}{}{0ex}{}{d+1}{i}\right)B_i$$
系数可以 $\mathrm{O}(d^2)$ 预处理。

则
$$\begin{array}{rl}{f}_d(n)& =\sum _{k\mid n}\mu (k)k^d{S}_d\left(\frac{n}{k}+1\right)\\ & =\sum _{k\mid n}\mu (k)k^d\sum _{i=0}^d{f}_i{\left(\frac{n}{k}+1\right)}^{d+1-i}\end{array}$$
你发现这个 ${\left(\frac{n}{k}+1\right)}^{m+1-i}$ 没法处理，但如果是 ${\left(\frac{n}{k}\right)}^{m+1-i}$ 就很好了——它可以和前面的 $k^d$ 相乘。

退回到上一步
$$\begin{array}{rl}{f}_d(n)& =\sum _{k\mid n}\mu (k)k^d\sum _{i=1}^{\frac{n}{k}}{i}^d\\ & =\sum _{k\mid n}\mu (k)k^d\sum _{i=1}^{\frac{n}{k}-1}{i}^d+\sum _{k\mid n}\mu (k)k^d{\left(\frac{n}{k}\right)}^d\\ & =\sum _{k\mid n}\mu (k)k^d\sum _{i=1}^{\frac{n}{k}-1}{i}^d+n^d\sum _{k\mid n}\mu (k)\\ & =\sum _{k\mid n}\mu (k)k^d\sum _{i=1}^{\frac{n}{k}-1}{i}^d+n^d[n=1]\end{array}$$
观察数据范围

$1\le \omega \le 10^3,2\le p_i\le 10^9,1\le {\alpha }_i\le 10^9$。

说明 $n\ne 1$。

则
$$\begin{array}{rl}{f}_d(n)& =\sum _{k\mid n}\mu (k)k^d\sum _{i=1}^{\frac{n}{k}-1}{i}^d\\ & =\sum _{k\mid n}\mu (k)k^d{S}_d\left(\frac{n}{k}\right)\\ & =\sum _{k\mid n}\mu (k)k^d\sum _{i=0}^d{f}_i{\left(\frac{n}{k}\right)}^{d+1-i}\\ & =\sum _{k\mid n}\mu (k)\sum _{i=0}^d{f}_i{n}^{d+1-i}{k}^{i-1}\\ & =\sum _{i=0}^d{f}_i{n}^{d+1-i}\sum _{k\mid n}\mu (k)k^{i-1}\end{array}$$
你会发现
$$\sum _{k\mid n}\mu (k)k^{i-1}=[(\mu \cdot {\mathrm{Id}}_{i-1})\ast 1](n)$$
所以这是一个积性函数，但是也不可能线性筛啊。

我们设
$$\sum _{k\mid n}\mu (k)k^{i-1}=g_i(n)$$
因为 $g_i(n)$ 是积性函数，所以可以按质因数拆开：
$$g_i(n)=\prod _{k=1}^{\omega }{g}_i(p_k^{{\alpha }_k})$$
则只用考虑质因数幂的情况。

你会发现
$$g_i(p^{\alpha })=\sum _{k\mid p^{\alpha }}\mu (k)k^{i-1}$$
其中只有 $k=1$ 和 $k=p$ 时 $\mu (k)$ 才有贡献。

而
$$\{\begin{array}{l}\mu (1)1^{i-1}=1\\ \mu (p)p^{i-1}=-p^{i-1}\end{array}$$
所以
$$g_i(p^{\alpha })=1-p^{i-1}$$
再贴一遍式子
$$f_d(n)=\sum _{i=0}^d{f}_i{n}^{d+1-i}\sum _{k\mid n}\mu (k)k^{i-1}$$
计算
$$g_i(n)=\prod _{k=1}^{\omega }1-p_k^{i-1}$$
是 $\mathrm{O}(\omega \log i)$ 的，求个和就是 $\mathrm{O}(\omega d\log d)$，但其实可以把 $p^{i-1}$ 往上递推做到 $\mathrm{O}(d\omega )$。

预处理伯努利数 $f_i$ 是 $\mathrm{O}(d^2)$ 的，而 $n^{d+1-i}$ 是 $\mathrm{O}(d\log d)$ 的，远不及前面的 $\mathrm{O}(d\omega )$。

综上，时间复杂度为 $\mathrm{O}(d^2+d\omega )$。

Code

注意预处理要处理到 $d+1$。

//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

const int MOD = 1e9 + 7;
const int MAXD = 105;

int add(int a, int b) {return (a + b) % MOD;}
int sub(int a, int b) {return (a - b + MOD) % MOD;}
int mul(int a, int b) {return (ll)a * b % MOD;}
int qpow(int a, int b) {int base = a, ans = 1; while (b) {if (b & 1) ans = mul(ans, base); base = mul(base, base); b >>= 1;} return ans;}
int GetInv(int a) {return qpow(a, MOD - 2);}

struct Bernoulli
{
	int fac[MAXD], fac_inv[MAXD], inv[MAXD];
	
	void pre(int d)
	{
		fac[0] = fac_inv[0] = inv[1] = 1;
		for (int i = 1; i <= d; i++)
		{
			fac[i] = mul(fac[i - 1], i);
		}
		fac_inv[d] = GetInv(fac[d]);
		for (int i = d - 1; i >= 1; i--)
		{
			fac_inv[i] = mul(fac_inv[i + 1], i + 1);
			inv[i + 1] = mul(fac_inv[i + 1], fac[i]);
		}
	}
	
	int C(int n, int m)
	{
		return mul(mul(fac[n], fac_inv[n - m]), fac_inv[m]);
	}
	
	int B[MAXD], f[MAXD];
	
	void cal(int d)
	{
		B[0] = 1;
		f[0] = inv[d + 1];
		for (int m = 1; m <= d; m++)
		{
			int sum = 0;
			for (int k = 0; k < m; k++)
			{
				sum = add(sum, mul(C(m + 1, k), B[k]));
			}
			B[m] = mul(MOD - sum, inv[m + 1]);
			f[m] = mul(mul(inv[d + 1], C(d + 1, m)), B[m]);
		}
	}
}Ber;

const int MAXW = 1005;

int p[MAXW], a[MAXW], xsy[MAXW];

int main()
{
	int d, w;
	scanf("%d%d", &d, &w);
	Ber.pre(d + 1), Ber.cal(d);
	int n = 1;
	for (int i = 1; i <= w; i++)
	{
		scanf("%d%d", p + i, a + i);
		n = mul(n, qpow(p[i], a[i]));
		xsy[i] = GetInv(p[i]);
	}
	int ans = 0;
	for (int i = 0; i <= d; i++)
	{
		int pro = 1;
		for (int j = 1; j <= w; j++)
		{
			pro = mul(pro, sub(1, xsy[j]));
			xsy[j] = mul(xsy[j], p[j]);
		}
		ans = add(ans, mul(mul(Ber.f[i], qpow(n, d + 1 - i)), pro));
	}
	printf("%d\n", ans);
	return 0;
}