bzoj 3601 一个人的数论

LINK：一个人的数论

这道题 是到好题。和伯努利数有关 但是我没学过。。

不难 把式子化简成\(\sum_{x|n}\mu(x)\cdot \sum_{i=1}^{\frac{n}{x}}(xi)^d\)

可以发现n巨大无比 我们除了能靠人类智慧拿一些分数之外就没办法了。

但是根据伯努利数 对于\(i^d\)求和是有一些办法的。

存在\(\sum_{i=1}^{n}i^d=\sum_{i=1}^{d+1}v_i n^i\) 其中\(v_i\)是常数。

那么原式还是可以化简的。

\(\sum_{i=1}^{d+1}v_i\sum_{x|n}\mu(x)\cdot x^d\cdot (\frac{n}{x})^i\)

可以发现 后面的东西是积性函数 我们可以对于每个p都求出来对应的答案最后再乘起来。

那么对于质因子\(p^a\)有\(\sum_{x|p^a}\mu(x)\cdot x^d\cdot (\frac{p^a}{x})^i\)

当\(x==1\)和\(x==p\)时才有值 所以总式=\(p^{ai}(1-p^{d-i})\)

考虑一下\(v_i\)怎么求。考虑使用伯努利数来求 就麻烦了 况且我也不会。

不过由于d只有100 可以直接列方程求解。上GAUSS即可。

考虑 一共有d+1个未知数 所以要列d+1个方程 等式的右边我们从1取到d+1即可。

从本题可以看出 题目中所给的约束条件都是有意义的 而不是无的放矢。

注意d==0时要保证自己的方程没有问题。

const ll MAXN=1010;
//$\sum_{i=1}^{n}i^d=\sum_{i=1}{d+1}v_i n^i$
ll a[110][110],f[MAXN],p[MAXN],v[MAXN];
ll n,d,ans;
inline ll ksm(ll b,ll p)
{
	ll cnt=1;
	while(p)
	{
		if(p&1)cnt=cnt*b%mod;
		b=b*b%mod;p=p>>1;
	}
	return cnt;
}
inline void GAUSS(ll n)
{
	rep(1,n,i)
	{
		ll p;
		rep(i,n,j)if(a[j][i]){p=j;break;}
		if(p!=i){rep(1,n,j)swap(a[i][j],a[p][j]);swap(f[i],f[p]);}
		rep(1,n,j)
		{
			if(i==j)continue;
			ll d=a[j][i]*ksm(a[i][i],mod-2)%mod;
			rep(1,n,k)a[j][k]=(a[j][k]-d*a[i][k])%mod;
			f[j]=(f[j]-f[i]*d)%mod;
		}
	}
	rep(1,n,i)
	{
		ll d=ksm(a[i][i],mod-2);
		f[i]=f[i]*d%mod;
	}
}
int main()
{
	freopen("1.in","r",stdin);
	get(d);get(n);
	rep(1,d+1,i)
	{
		ll w=1;a[i][0]=1;
		rep(1,d+1,j)
		{
			w=w*i%mod;
			a[i][j]=w;
		}
		f[i]=f[i-1]+a[i][d];
	}
	GAUSS(d+1);
	rep(1,n,i)get(p[i]),get(v[i]);
	if(n==1&&p[1]==1){puts("1");return 0;}
	rep(1,d+1,i)
	{
		ll cnt=1;
		rep(1,n,j)
			cnt=cnt*ksm(p[j],(v[j]*i)%(mod-1))%mod*(1-ksm(p[j],(mod-1+d-i)%(mod-1)))%mod;
		ans=(ans+f[i]*cnt%mod)%mod;
	}
	putl((ans+mod)%mod);
	return 0;
}