[BZOJ3601]一个人的数论-题解

【题目地址】

题意简述

给你$n,d$，求下面式子在$\mathrm{m}\mathrm{o}\mathrm{d}10^9+7$意义下的值。

$$\sum _{i=1}^n[gcd(i,n)=1]i^d$$

其中$n$特别大，所以我们给你一个$w$，然后给你两个数组$p_i,c_i$，其中：

$$n=\prod _{i=1}^w{p}_i^{c_i}$$

$d\le 100,w\le 1000,p_i,c_i\le 10^9$

---

我们仍旧先推一波式子：

$$\begin{array}{rl}ans=& \sum _{i=1}^n[gcd(i,n)=1]i^d\\ =& \sum _{i=1}^n{i}^d\sum _{j|i,j|n}\mu (j)\\ =& \sum _{j=1}^n\mu (j)\sum _{j|i}^n{i}^d\\ =& \sum _{j=1}^n\mu (j)\sum _{i=1}^{\lfloor \frac{n}{j}\rfloor }(ij{)}^d\\ =& \sum _{j=1}^n\mu (j)j^d\sum _{i=1}^{\lfloor \frac{n}{j}\rfloor }{i}^d\end{array}$$

后面部分就是自然数幂的前缀和，拉格朗日差值就可以了，在$O(dlogd)$时间内求出，而前面是可以杜教筛的：

$$\begin{array}{rl}\sum _{j|n}\mu (j)j^d=& (\mu \cdot i{d}^d)\ast 1\end{array}$$

我们令$f=(\mu \cdot i{d}^d)\ast 1,g=i{d}^d$，然后$f\ast g$就等于：

$$\sum _{j|n}\mu (j)j^d(\frac{n}{j}{)}^d=n^d\sum _{j|n}\mu (j)=n^d[n=1]=[n=1]=ϵ$$

所以$f=(\mu \cdot i{d}^d)\ast 1,g=i{d}^d,h=ϵ$，然后套入杜教筛公式即可：

$$S(n)=1-\sum _{i=2}^n{i}^{d}S(\lfloor \frac{n}{i}\rfloor )$$

然后就可以在$O(n^{\frac{2}{3}}dlogd)$的时间内解决了。

---

等等，$n$好像有$10^{10000...}$左右，-_-||。

但是不是就这样束手就擒，打暴力去了呢？

不，我们离正解还差一步了。

---

但是，我们根据拉格朗日差值的经验可以知道，${\sum }_{i=1}^n{i}^d$为一个$d+1$次多项式，那么我们令$g(x)={\sum }_{i=1}^x{i}^d$，同时我们也可以知道$g(x)={\sum }_{i=1}^{d+1}{a}_i{x}^i$，假设我们求出了所有$a_i$，我们将其带入式子，看看有没有什么奇妙的反应：

$$\begin{array}{rl}ans=& \sum _{j=1}^n\mu (j)j^d\sum _{i=1}^{\lfloor \frac{n}{j}\rfloor }{i}^d\\ =& \sum _{j=1}^n\mu (j)j^d\sum _{i=1}^{d+1}{a}_i{\left(\lfloor \frac{n}{j}\rfloor \right)}^i\\ =& \sum _{i=1}^{d+1}{a}_i\sum _{j|n}\mu (j)j^d{\left(\lfloor \frac{n}{j}\rfloor \right)}^i\end{array}$$

现在我们主要是要解决后半部分的快速求法，我们令$f(i)={\sum }_{j|n}\mu (j)j^d{\left(\lfloor \frac{n}{j}\rfloor \right)}^i$，那么我们知道$f=(\mu \cdot i{d}^d)\ast i{d}^i$，那么根据$\mu ,i{d}^k$都为积性函数，而积性函数的狄利克雷卷积都为积性函数，所以我们可以知道$f$也是一个积性函数。

这里，我们看，题目都将$n$的质因数帮我们分解好了凯森(ノ￣▽￣)，所以我们将质因数分开考虑，因为$f(n)=\prod f(p_i^{c_i})$，所以此时，我们只需知道$f(p^c)$如何计算：

我们看原式
$$\sum _{j|p^c}\mu (j)j^d{\left(\lfloor \frac{n}{j}\rfloor \right)}^i=\sum _{k=0}^c\mu (p^k)(p^k{)}^d{\left(\lfloor \frac{n}{p^k}\rfloor \right)}^i$$

而根据$\mu$的定义，我们知道当$k>1$时的时候$\mu (p^k)=0$，所以$k=0,1$，那么$f(p^c)$就等于：

$$(p^c{)}^i-p^d(p^{c-1}{)}^i=p^{ci}-p^{ci-i+d}$$

那么枚举$a_i$，然后枚举$p_j$，然后算就用快速幂，复杂度就是$O(dwlog(10^9+6))$（因为根据欧拉定理和费马小定理，我们可以将指数$\mathrm{m}\mathrm{o}\mathrm{d}10^9+6$）。

好啦，就解决啦！ヾ(@^▽^@)ノ

---

等等辣QAQ，$a_i$好像还没说怎么求-_-||。

其实，我们知道它是一个$d+1$次多项式的系数，那么高斯消元就好啦，我们列出$g(1)\sim g(d+1)$的方程即可。

前面回顾：$g(x)={\sum }_{i=1}^{d+1}{a}_i{x}^i$

这个方程组的初始化，我们每次枚举$x$，然后用快速幂计算出当前的$g(x)$的值，然后枚举$a_i$，依次计算出$x^i$作为系数即可，然后就是高斯消元板子啦~~

---

丑陋的加调了半天的代码_(¦3」∠)_：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=1e3+10,M=1e2+10;
const ll Mod=1e9+7;
int d,w;
ll fpow(ll a,ll b){
	ll res=1;
	for(;b;b>>=1,a=(a*a)%Mod)
		if(b&1)res=(res*a)%Mod;
	return res;
}
ll A[M][M],X[M],P[N],C[N];
void Guass(){
	int D=d+1;
	for(ll i=0,sum=0,x;i<=D;i++){
		x=1;sum=(sum+fpow(i,d))%Mod;
		A[i][D+1]=sum;
		for(ll j=0;j<=D;j++){
			A[i][j]=x;x=(x*i)%Mod;	
		}
	}
	for(int now=0,pos;now<=D;now++){
		pos=now;
		if(!A[pos][now])
			for(int i=now+1;i<=D;i++)
				if(A[i][now]){pos=i;break;}
		if(pos!=now)swap(A[now],A[pos]);
		ll inv=fpow(A[now][now],Mod-2);
		for(int i=now+1;i<=D;i++){
			ll tt=(A[i][now]*inv)%Mod;
			for(int j=0;j<=D+1;j++){
				A[i][j]=(A[i][j]-A[now][j]*tt%Mod)%Mod;
				if(A[i][j]<0)A[i][j]+=Mod;
			}
		}
	}
	for(int i=D;i>=0;i--){
		for(int j=D;j>=i+1;j--){
			A[i][D+1]=(A[i][D+1]-A[i][j]*X[j]%Mod)%Mod;
			if(A[i][D+1]<0)A[i][D+1]+=Mod;
		}
		X[i]=A[i][D+1]*fpow(A[i][i],Mod-2)%Mod;
	}
}
ll GetH(ll a,ll b){
	ll p1=C[b]*a%(Mod-1);
	ll p2=(p1+d-a)%(Mod-1);if(p2<0)p2+=(Mod-1);
	ll ans=(fpow(P[b],p1)-fpow(P[b],p2))%Mod;
	if(ans<0)ans+=Mod;
	return ans;
}
ll ans;
int main(){
	scanf("%d%d",&d,&w);
	Guass();
	for(int i=1;i<=w;i++)scanf("%lld%lld",&P[i],&C[i]);
	for(int i=0,D=d+1;i<=D;i++){
		ll xx=X[i];
		for(int j=1;j<=w;j++)xx=(xx*GetH(i,j))%Mod;
		ans=(ans+xx)%Mod;
	}
	printf("%lld\n",ans);
	return 0;
}

---

End

元旦节只放一天，1月1日还要在学校度过，惨兮兮QWQ