【LOJ6207】米缇（莫比乌斯反演）（杜教筛）（拉格朗日插值）

传送门

题解：

杜教筛里面long long 没有取模炸了一堆。。。

首先推式子：

A n s = ∑ i = 1 n ∑ j = 1 n σ k ( i j ) = ∑ i = 1 n ∑ j = 1 n ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = 1 ] ( i y x ) k = ∑ d = 1 n μ ( d ) d k ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ ∑ x ∣ i ∑ y ∣ j ( i y x ) k = ∑ d = 1 n μ ( d ) d k ∑ i = 1 ⌊ n d ⌋ ∑ x ∣ i ( i x ) k ∑ j = 1 ⌊ n d ⌋ ∑ y ∣ j y k = ∑ d = 1 n μ ( d ) d k ( ∑ i = 1 ⌊ n d ⌋ i k ⌊ n i d ⌋ ) 2 \begin{aligned} Ans&=&&\sum_{i=1}^{n}\sum_{j=1}^n\sigma_k(ij)\\ &=&&\sum_{i=1}^n\sum_{j=1}^n\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)=1](\frac{iy}{x})^k\\ &=&&\sum_{d=1}^n\mu(d)d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{x\mid i}\sum_{y\mid j}(\frac{iy}{x})^k\\ &=&&\sum_{d=1}^n\mu(d)d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{x\mid i}(\frac{i}{x})^k\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y\mid j}y^k\\ &=&&\sum_{d=1}^n\mu(d)d^k(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i^k\lfloor\frac{n}{id}\rfloor)^2 \end{aligned} Ans​=====​​i=1∑n​j=1∑n​σk​(ij)i=1∑n​j=1∑n​x∣i∑​y∣j∑​[gcd(x,y)=1](xiy​)kd=1∑n​μ(d)dki=1∑⌊dn​⌋​j=1∑⌊dn​⌋​x∣i∑​y∣j∑​(xiy​)kd=1∑n​μ(d)dki=1∑⌊dn​⌋​x∣i∑​(xi​)kj=1∑⌊dn​⌋​y∣j∑​ykd=1∑n​μ(d)dk(i=1∑⌊dn​⌋​ik⌊idn​⌋)2​

对于 μ ⋅ I d k \mu\cdot Id_k μ⋅Idk​，我们利用 I d k Id_k Idk​可以直接构造杜教筛

对于后面的 S σ k S_{\sigma_k} Sσk​​，可以直接筛一部分后暴力做。

对于 I d k Id_k Idk​，直接用拉格朗日插值法算一算就行了。

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;

cs int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}
inline void Inc(int &a,int b){(a+=b)>=mod&&(a-=mod);}
inline void Dec(int &a,int b){(a-=b)<0&&(a+=mod);}
inline int MD(ll a){return a<0?(a<-mod?(a%mod+mod)%mod:a+mod):(a>=mod?a%mod:a);}

cs int N=1e7,K=7e3,SIZE=N+5;

ll n;int k,lim;

bool mark[SIZE];
int dk[SIZE],mu[SIZE],dm[SIZE],pw[SIZE];
std::vector<int> prime;
int fac[K+5],ifac[K+5];

inline void init(int len){lim=len;
	mu[1]=pw[1]=dm[1]=dk[1]=1;
	for(int re i=2;i<=len;++i){
		if(!mark[i]){
			prime.push_back(i);
			dm[i]=1;
			pw[i]=power(i,k);
			dk[i]=add(pw[i],1);
			mu[i]=mod-1;
		}
		for(int re j:prime){
			if(i*j>len)break;
			mark[i*j]=true;
			pw[i*j]=mul(pw[i],pw[j]);
			if(i%j){
				dm[i*j]=i;
				mu[i*j]=dec(0,mu[i]);
				dk[i*j]=mul(dk[i],dk[j]);
			}
			else {
				dm[i*j]=dm[i];
				dk[i*j]=dm[i]==1?add(1,mul(dk[i],pw[j])):mul(dk[i/dm[i]*j],dk[dm[i]]);
				break;
			}
		}
	}
	for(int re i=2;i<=len;++i){
		mu[i]=add(mu[i-1],mul(mu[i],pw[i]));
		Inc(dk[i],dk[i-1]);
		Inc(pw[i],pw[i-1]);
	}
	fac[0]=fac[1]=ifac[0]=1;
	for(int re i=2;i<=k+1;++i)fac[i]=mul(fac[i-1],i);
	ifac[k+1]=power(fac[k+1],mod-2);
	for(int re i=k;i;--i)ifac[i]=mul(ifac[i+1],i+1);
}

struct Map{
	static cs int magic=18985739;
	ll key[magic];
	int val[magic],tmp;
	Map(){memset(key,-1,sizeof key);}
	inline int &operator[](ll k){
		int h=k%magic;
		while(key[h]!=-1&&key[h]!=k)h=h+1==magic?0:h+1;
		if(key[h]==-1){key[h]=k;val[h]=0;}
		return val[h];
	}
	inline cs int &operator[](ll k)cs{
		int h=k%magic;
		while(key[h]!=-1&&key[h]!=k)h=h+1==magic?0:h+1;
		return key[h]==-1?tmp:val[h];
	}
	inline bool find(ll k)cs{
		int h=k%magic;
		while(key[h]!=-1&&key[h]!=k)h=h+1==magic?0:h+1;
		return key[h]==k;
	}
}sum_k;

int pre[K+5],suf[K+5];
inline int Sum_k(ll n){
	if(n<=lim)return pw[n];
	if(sum_k.find(n))return sum_k[n];
	pre[0]=suf[k+3]=1;int t=n%mod;
	for(int re i=1;i<=k+2;++i)pre[i]=mul(pre[i-1],dec(t,i));
	for(int re i=k+2;i>=1;--i)suf[i]=mul(suf[i+1],dec(t,i));
	int ans=0;
	for(int re i=1;i<=k+2;++i){
		int coef=mul(mul(pre[i-1],suf[i+1]),mul(ifac[i-1],ifac[k+2-i]));
		(k^i)&1?Dec(ans,mul(coef,pw[i])):Inc(ans,mul(coef,pw[i]));
	}
	return sum_k[n]=ans;
}

inline int Sum_dk(ll n){
	if(n<=lim)return dk[n];
	int ans=0;
	for(ll i=1,j,last=0,now;i<=n;i=j+1){
		now=Sum_k(j=n/(n/i));
		Inc(ans,mul(dec(now,last),MD(n/i)));
		last=now;
	}
	return ans;
}

bool vis[100005];
int f[100005];
inline int Sum_mu(ll n){
	if(n<=lim)return mu[n];
	if(vis[::n/n])return f[::n/n];
	vis[::n/n]=true;
	int ans=1;
	for(ll re i=2,j,last=1,now;i<=n;i=j+1){
		now=Sum_k(j=n/(n/i));
		Dec(ans,mul(dec(now,last),Sum_mu(n/i)));
		last=now;
	}
	return f[::n/n]=ans;
}

signed main(){
//	freopen("miti.in","r",stdin);
	std::cin>>n>>k;
	init(std::min(N*1.,std::max(k*1.+1,pow(n,2./3.))));
	int ans=0;
	for(ll re i=1,j,last=0,now;i<=n;i=j+1){
		j=n/(n/i);
		int sum=Sum_dk(n/i);
		now=Sum_mu(j);
		Inc(ans,mul(dec(now,last),mul(sum,sum)));
		last=now;
	}
	cout<<ans<<"\n";
	return 0;
}