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

一、题目

点此看题

二、解法

你会发现好像没有什么巧妙的算法，不如我们直接把答案形式化地写出来：

\[\sum_{i=1}^n[\gcd(i,n)=1]\space i^k \]

然后看到了熟悉的 \(\gcd\) 结构，我们直接反演：

\[\sum_{i=1}^ni^k\sum_{x|(i,n)}\mu(x) \]

然后我们枚举 \(x\) ：

\[\sum_{x|n}\mu(x)x^k\sum_{i=1}^{n/x} i^k \]

你会发现还是很复杂，现在的契机在于后面的 \(\sum_{i=1}^{n/x}i^k\) 是有很多方法的，但是我们的标准只有一个：把 \(k\) 融入到枚举中所以我们就能得到关于 \(k\) 的复杂度 ，我们可以试一试用第二类斯特林数展开它（但是我试过不行），还有一种方法，因为他是一个 \(k+1\) 次的关于 \(n/x\) 的多项式，我们假设已经计算了他的系数 \(f_i\) ，这样改写算式：

\[\sum_{x|n}\mu(x)x^k\sum_{i=1}^kf_i\times(\frac{n}{x})^i \]

然后我们先枚举 \(i\) ，这也是我们希望看到的：

\[\sum_{i=1}^kf_i\sum_{x|n}\mu(x)x^k(\frac{n}{x})^i \]

\[\sum_{i=1}^kf_in^i\sum_{x|n}\mu(x)x^{k-i} \]

现在要好好观察后面一部分了，他是我们的复杂度瓶颈，这道题很特殊的地方是题目给了我们一个乘积式来表示 \(n\) ，后面那一部分又一定是一个积性函数（因为 \(\mu\) 和 \(x^{k-i}\) 都是可以拆分的），所以可以对每个质数分开考虑。

对于质数 \(p_i^{a_i}\) ，当 \(x=1\) 时答案是 \(1\) ，当 \(x=p_i\) 时答案是 \(-p_i^{k-i}\) ，对于其他情况 \(\mu=0\) 所以不用考虑。

---

现在问题变成了算系数 \(f_i\) ，你可以用高斯消元爆操，但是有一种更优美的方法：伯努利数。伯努利数就是用来处理这种幂求和的问题，学会这种方法你只需要记住两个公式，\(B_i\) 表示伯努利的定义式是这样的：

\[\sum_{i=0}^n C(n+1,i)B_i=[n=0] \]

那么我们 \(O(k^2)\) 就可以求 \(B_i\) ，怎么求 \(f_i\) 呢？首先要知道伯努利公式：

\[\sum_{i=1}^ni^k=\frac{1}{k+1}\sum_{i=0}^kC(k+1,i)B_in^{k+1-i} \]

所以我们就知道了 \(f_{k+1-i}=\frac{C(k+1,i)B_i}{k+1}\) ，证明我暂时不懂，以后再填坑吧。

#include <cstdio>
const int M = 1005;
const int MOD = 1e9+7;
#define int long long
int read()
{
	int x=0,f=1;char c;
	while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
	while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
int k,w,ans,p[M],a[M],B[M],f[M],fac[M],inv[M];
int C(int n,int m)
{
	if(n<m) return 0;
	return fac[n]*inv[m]%MOD*inv[n-m]%MOD;
}
int qkpow(int a,int b)
{
	int r=1;
	while(b>0)
	{
		if(b&1) r=r*a%MOD;
		a=a*a%MOD;
		b>>=1;
	}
	return r;
}
void init()
{
	fac[0]=inv[0]=inv[1]=1;
	for(int i=1;i<=k+1;i++) fac[i]=fac[i-1]*i%MOD;
	for(int i=2;i<=k+1;i++) inv[i]=inv[MOD%i]*(MOD-MOD/i)%MOD;
	for(int i=2;i<=k+1;i++) inv[i]=inv[i-1]*inv[i]%MOD;
	B[0]=1;
	for(int i=1;i<=k;i++)
	{
		int sum=0;
		for(int j=0;j<i;j++)
			sum=(sum-C(i+1,j)*B[j])%MOD;
		B[i]=sum*qkpow(i+1,MOD-2)%MOD;
	}
	for(int i=0;i<=k;i++)
		f[k+1-i]=C(k+1,i)*B[i]%MOD*qkpow(k+1,MOD-2)%MOD;
}
signed main()
{
	k=read();w=read();
	for(int i=1;i<=w;i++)
	{
		p[i]=read();a[i]=read();
	}
	init();
	for(int i=1;i<=k+1;i++)
	{
		int res=f[i];
		for(int j=1;j<=w;j++)
		{
			res=res*qkpow(p[j],a[j]*i)%MOD;//这里我懒了，多了个log 
			res=res*(1-qkpow(p[j],k-i+MOD-1))%MOD;
		}
		ans=(ans+res)%MOD;
	}
	printf("%lld\n",(ans+MOD)%MOD);
}

标签：frac,队互测,数论,sum,mu,int,2014,我们,伯努利
来源： https://www.cnblogs.com/C202044zxy/p/14264615.html