bzoj 3601: 一个人的数论（莫比乌斯反演 + 高斯消元 + 积性函数）

---

题目要求的是 $f(n,d)=\sum _{i=1}^n[gcd(i,n)=1]i^d$

设 $$h(x)=\sum _{i=1}^n[gcd(i,n)=x]i^d$$$$H(x)=\sum _{x|p}h(p)$$

（从贡献角度）等价于 $H(x)=\sum _{i=1}^n[x|gcd(i,n)]i^d$
根据莫比乌斯反演，$h(x)=\sum _{x|p}H(p)\ast \mu (\frac{p}{x})$，得到
$$h(1)=\sum _{p=1}^n\mu (p)\ast \sum _{i=1}^n[p|gcd(i,n)]i^d$$$$h(1)=\sum _{p|n}\mu (p)\ast \sum _{i=1}^{\frac{n}{p}}(ip{)}^d=\sum _{p|n}\mu (p)\ast p^d\sum _{i=1}^{\frac{n}{p}}{i}^d$$
到这一步似乎无法做下去了，因为 n 特别大。此时只能根据题解大胆猜测，数据这么大还给的是质因子分解的形式，可能要用积性函数来做，因为积性函数满足 $f(x\ast y)=f(x)\ast f(y)\ast [gcd(x,y)=1]$，可以计算每一项质因子的贡献最后乘起来得到答案。

观察发现这个式子是一个接近卷积的形式，前两个函数都是积性函数，根据经验可以得知后面一个求和式是一个以 $\frac{n}{p}$ 为自变量的 $d+1$ 次多项式，但这个多项式并不是一个积性函数。

注意到要将整个式子凑成一个积性函数，右边那个和式必须是一个以 $\frac{n}{p}$ 为变量的积性函数，这样整个式子可以凑成一个 Dirichlet 卷积，进一步根据题解大胆尝试，将右边$d+1$ 以多项式的形式传入参数 $\frac{n}{p}$ 代入得到：$$h(1)=\sum _{p|n}\mu (p)\ast p^d\sum _{i=0}^{d+1}{b}_i\ast (\frac{n}{p}{)}^i$$改变枚举项，先枚举 $i$，得到$$h(1)=\sum _{i=0}^{d+1}{b}_i\sum _{p|n}\mu (p)\ast p^d(\frac{n}{p}{)}^i$$

这个时候右边那部分就是积性函数了，并且由于 $\mu (x)$ 函数的特性，每一项 $p_i^{a_i}$ 只有 1 和 pi 有贡献，只要计算出这个多项式的所有系数 bi 就可以在规定时间内求解。

对于多项式求系数 bi 有两种方法：拉格朗日插值多点插值，高斯消元

拉格朗日多点插值求解多项式系数还不会，且这里 $d\le 100$，采用高斯消元进行求解

---

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
typedef long long ll;
const int maxn = 1e3 + 10;
int k,w;
int p[maxn],a[maxn],b[maxn];
int x[maxn][maxn],y[maxn];
inline int read() {
	int x=0,f=1;char ch;
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}

inline int add(int x, int y) {
  	x += y;
  	if (x >= mod) {
    	x -= mod;
  	}
  	return x;
}

inline int sub(int x, int y) {
  	x -= y;
  	if (x < 0) {
    	x += mod;
  	}
  	return x;
}

inline int mul(int x, int y) {
  	return (long long) x * y % mod;
}
int fpow(int a,int b) {
	int r = 1;
	while (b) {
		if (b & 1) r = mul(r,a);
		a = mul(a,a);
		b >>= 1;
	}
	return r;
}
void Gauss(int x[maxn][maxn],int y[maxn],int n) {
	for (int i = 0; i <= n; i++) {
		for (int j = 0; j < i; j++) {
			int tp = mul(fpow(x[i][i],mod - 2),x[j][i]);
			for (int k = 0; k <= n; k++) {
				x[j][k] = sub(x[j][k],mul(tp,x[i][k]));
			}
			y[j] = sub(y[j],mul(tp,y[i]));
		}
		for (int j = i + 1; j <= n; j++) {
			int tp = mul(fpow(x[i][i],mod - 2),x[j][i]);
			for (int k = 0; k <= n; k++) {
				x[j][k] = sub(x[j][k],mul(tp,x[i][k]));
			}
			y[j] = sub(y[j],mul(tp,y[i]));
		}
	}
}
int main() {
	k = read(); w = read();
	for (int i = 1; i <= w; i++)
		p[i] = read(), a[i] = read();
	y[0] = 0;
	for (int i = 0; i <= k + 1; i++) {
		for (int j = 0; j <= k + 1; j++)
			x[i][j] = fpow(i,j);
		if (i > 0) y[i] = add(y[i - 1],fpow(i,k));
	}
	Gauss(x,y,k + 1);
	for (int i = 0; i <= k + 1; i++)
		b[i] = mul(y[i],fpow(x[i][i],mod - 2));
	int sum = 0, ans = 1;
	for (int i = 0; i <= k + 1; i++) {
		ans = 1; 
		for (int j = 1; j <= w; j++) {
			int res = 0;
			int t1 = fpow(p[j],a[j]), t2 = fpow(p[j],a[j] - 1);
			res = add(res,fpow(t1,i));
			res = sub(res,mul(fpow(t2,i),fpow(p[j],k)));
			ans = mul(ans,res);
		}	
		sum = add(sum,mul(b[i],ans));
	}
	printf("%d\n",sum);
	return 0;
}