bzoj3601-一个人的数论

题目

定义：

\[\begin{aligned} f_d(n)=\sum _{gcd(x,n)=1} x^d \end{aligned} \]

给定\(d\)和\(n\) (\(n\)以质因数分解的形式给出)，求\(f_d(n)\ (mod \ 10^9+7)\)。

分析

本来是学莫比乌斯反演时候的题，拖到了高斯消元来做。
看到一个互质，即\(gcd(x,n)=1\)，马上就能想到莫比乌斯函数的那个性质：

\[令 \ f(n)=\sum _{d|n} \mu(d) \\ 则有 \ f(n)= \left\{\begin{aligned} 1，n=1\\0，n>1\end{aligned} \right. \]

下面我们证明这个结论：
若\(n=1\)，结论显然。
若\(n>1\)，那么我们设

\[\begin{aligned} s=\prod _{i=1} ^k p_i^{a_i}，k|n \end{aligned} \]

当任意一个\(a\_i>1\)时，\(\mu(s)=0\)
所以上面的\(f(n)\)等价于\(f(s)\)，其中\(s=\prod\_{i=1}^k p\_i\)
根据莫比乌斯函数的定义，当\(k\)为奇数时，\(f(s)=-1\)，\(k\)为偶数时，\(f(s)=1\).
所以\(f(n)\)可以写成：

\[\begin{aligned} f(n)=C_k^0-C_k^1+C_k^2....=\sum _{i=0}^k (-1)^i*C_k^i \end{aligned} \]

所以我们现在要证明的就是：

\[\begin{aligned} \sum _{i=0}^k (-1)^i*C_k^i=0 \end{aligned} \]

注意到二项式定理：

\[\begin{aligned} (x+y)^n=\sum _{i=0}^k C_k^i*x^i*y^{k-i} \end{aligned} \]

与我们要证明的结论的相似性，我们令\(x=-1\)，\(y=1\)，即可得证。

这样，我们可以对原式进行变形：

\[\begin{aligned} f_d(n)&=\sum _{gcd(x,n)=1} x^d \\ &=\sum _{x=1}^n x^d[gcd(x,n)=1] \\ &=\sum _{x=1}^n x^d\sum _{c|gcd(x,n)}\mu(c) \\ &=\sum _{c|n} \sum _{x=1}^{\frac{n}{c}} (cx)^d\mu(c) \\ &=\sum _{c|n} \mu(c)c^d \sum _{x=1} ^{\frac{n}{c}} x^d \\ \end{aligned} \]

令：

\[\begin{aligned} h(n)=\sum _{i=1}^n i^d \end{aligned} \]

我们称\(h\)为自然数幂和。
则有：

\[\begin{aligned} f_d(n)=\sum _{c|n} \mu(c)c^d h(\frac{n}{c}) \\ \end{aligned} \]

这里有一个很奇妙的，并且正确的猜想：

\[\begin{aligned} h(n)=\sum _{i=1} ^{d+1}a_in^i \end{aligned} \]

也就是说，自然数幂和可以用\(n\)的多项式表示出来。
于是奇妙地用到了高斯消元，直接用\(1\)到\(d+1\)为例子，列出\(d+1\)个方程，强行把多项式的系数解出来。其实自然数幂和有通项公式，近期会有相关的文章。于是我们有：

\[\begin{aligned} f_d(n)&=\sum _{c|n} \mu(c)c^d \sum _{i=1} ^{d+1}a_i(\frac{n}{c})^i \\ &=\sum _{i=1}^{d+1} a_i \sum _{c|n} \mu(c)c^d(\frac{n}{c})^i \end{aligned} \]

令：

\[\begin{aligned} g_i(n)=\sum _{c|n} \mu(c)c^d(\frac{n}{c})^i \end{aligned} \]

那么有：

\[\begin{aligned} f_d(n)=\sum _{i=1}^{d+1} a_i *g_i(n) \end{aligned} \]

现在我们可以在\(O(d^3)\)的时间内求出所有的\(a_i\)，只要对每个\(i\)求出\(g_i(n)\)，我们就能\(O(d)\)算出答案了。
若令\(r(c)=\mu(c)c^d\)，那么有：

\[\begin{aligned} g_i(n)=\sum _{c|n} r(c)(\frac{n}{c})^i \end{aligned} \]

这是一个狄利克雷卷积(即形如 \(t(n)=\sum _{d|n}f(d)g(\frac{n}{d})\)的函数)！
狄利克雷卷积函数有一个性质，如果原来两个函数都是积性函数，那么新函数也是积性函数。证明如下：
如果两个函数不完全积性，那么需要\(gcd(x,y)=1\)，否则不需要。

\[\begin{aligned} t(x)t(y)&=\sum _{d|x}f(d)g(\frac{x}{d}) \sum _{e|y}f(e)g(\frac{y}{e}) \\ &=\sum _{de|xy}f(de)g(\frac{xy}{de}) \\ &=t(xy) \end{aligned} \]

所以我们可以通过\(n\)的质因数分解求出每个\(g_i(n)\)：

\[\begin{aligned} g_i(n)=\prod _{j=1}^w g_i(p_j^{q_j}) \end{aligned} \]

而对于只包含一个质数的\(g\)函数：

\[\begin{aligned} g_i(p^q) &=\sum _{j=0}^q \mu(p^j)p^{jd}(\frac{p^q}{p^j})i \\ &= \left\{\begin{aligned} p^{qi}，j=0\\-p^{qi+d-i}，j=1\\ 0, otherwise\end{aligned} \right. \end{aligned} \]

所以有：

\[\begin{aligned} g_i(p^q) =p^{qi}-p^{qi+d-i} \end{aligned} \]

直接计算即可。还是很简单的嘛！

代码

公式推错了四次～

#include<cstdio>
#include<algorithm>
#include<cctype>
using namespace std;
typedef long long giant;
giant read() {
	giant x=0,f=1;
	char c=getchar();
	for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*f;
}
const giant maxn=1e3+5;
const giant maxd=105;
const giant q=1e9+7;
giant p[maxn],s[maxn];
giant a[maxd][maxd],as[maxn];
giant mi(giant x,giant y) {
	giant ret=1;
	while (y) {
		if (y&1) (ret*=x)%=q;
		y>>=1,(x*=x)%=q;
	}
	return ret;
}
void elm(giant n) {
	for (giant i=1;i<n;++i) {
		if (!a[i][i]) {
			for (giant j=i+1;j<=n;++j) if (a[j][i]) {
				for (giant k=1;k<=n+1;++k) swap(a[j][k],a[i][k]);
				break;
			} 	
		}	
		giant inv=mi(a[i][i],q-2);
		for (giant j=i+1;j<=n;++j) if (a[j][i]) {
			giant tmp=(a[j][i]*inv)%q;
			for (giant k=1;k<=n+1;++k) a[j][k]=((a[j][k]-a[i][k]*tmp+q)%q+q)%q;
		}
	}
	for (giant i=n;i>1;--i) {
		if (!a[i][i]) {
			for (giant j=i-1;j>0;--j) if (a[j][i]) {
				for (giant k=1;k<=n+1;++k) swap(a[j][k],a[i][k]);
				break;
			}	
		}	
		giant inv=mi(a[i][i],q-2);
		for (giant j=i-1;j>0;--j) if (a[j][i]) {
			giant tmp=(a[j][i]*inv)%q;
			for (giant k=1;k<=n+1;++k) a[j][k]=((a[j][k]-a[i][k]*tmp+q)%q+q)%q;	
		}
	}
}
int main() {
	#ifndef ONLINE_JUDGE
		freopen("test.in","r",stdin);
		freopen("my.out","w",stdout);
	#endif
	giant d=read(),w=read();
	for (giant i=1;i<=w;++i) p[i]=read(),s[i]=read();
	giant sum=0;
	for (giant i=1;i<=d+1;++i) {
		(sum+=mi(i,d))%=q;
		giant tmp=1;
		for (giant j=1;j<=d+1;++j) (tmp*=i)%=q,a[i][j]=tmp;
		a[i][d+2]=sum;
	}
	elm(d+1); 
	for (giant i=1;i<=d+1;++i) as[i]=(a[i][d+2]*mi(a[i][i],q-2))%q;
	giant ans=0;
	for (giant i=1;i<=d+1;++i) {
		giant g=1;
		for (giant j=1;j<=w;++j) {
			giant tmp=(mi(p[j],s[j]*i)-mi(p[j],s[j]*i+d-i)+q)%q;
			(g*=tmp)%=q;
		}
		(ans+=(as[i]*g))%=q;
	}
	printf("%lld\n",ans);
}