【HDU6607】Easy Math Problem（杜教筛）（min_25筛）（拉格朗日插值）

传送门

---

题解：

推起来还挺简单的一道题。

我们要求的是：

$$\begin{array}{rl}Ans=& \sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j{)}^{k}lcm(i,j)[gcd(i,j)\in P]\\ =& \sum _{p\in P}\sum _{i=1}^n\sum _{j=1}^n{p}^{k-1}ij[gcd(i,j)=p]\\ =& \sum _{p\in P}{p}^{k+1}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }i\sum _{j=1}^{\lfloor \frac{n}{p}\rfloor }j[gcd(i,j)=1]\end{array}$$

这个推了不知道多少次的东西。。。我们知道：

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

这个玩意卷上一个 $i{d}^2$ 即可杜教筛。

于是还剩下的问题就是求出区间质数的 $k+1$ 次方之和，直接拉格朗日插值+min_25筛即可。

复杂度为 $O(\sqrt{n}k+n^{3/4}/\log n)$

---

代码：

#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,iv6=(mod+1)/6,iv2=(mod+1)/2;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(int a,int b){return a-b<0?a-b+mod:a-b;}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int po(int a,int b){int r=1;for(;b;b>>=1,Mul(a,a))if(b&1)Mul(r,a);return r;}
inline void ex_gcd(int a,int b,int &x,int &y){
	if(!b){x=1,y=0;return;}ex_gcd(b,a%b,y,x);y-=a/b*x;
}inline int Inv(int a){int x,y;ex_gcd(mod,a,y,x);return x+(x>>31&mod);}
inline int sqr(int x){return mul(x,x);}
inline int cub(int x){return mul(x,mul(x,x));}
inline int MD(ll x){return x>=mod?x%mod:x;}
inline int S1(ll n){n=MD(n);return MD(n*(n+1)>>1);}
inline int S2(ll n){n=MD(n);return mul(mul(n,n+1),mul(n<<1|1,iv6));}
inline int S3(ll n){n=MD(n);return sqr(MD(n*(n+1)>>1));}

cs int N=5e6+7;
cs int M=2e3+7;
cs int K=1e2+7;
cs int NN=1e5+7;

ll n;
int k,tim,nn;
int p[N],cp;
bool mrk[N];
int f[N],F[M];

void linear_sieves(){
	f[1]=1;
	for(int re i=2;i<N;++i){
		if(!mrk[i])
			p[++cp]=i,f[i]=mul(sqr(i),i-1);
		for(int re j=1;i*p[j]<N;++j){
			mrk[i*p[j]]=true;
			if(i%p[j]){
				f[i*p[j]]=mul(f[i],f[p[j]]);
			}else {
				f[i*p[j]]=mul(f[i],cub(p[j]));break;
			}
		}
	}
	for(int re i=2;i<N;++i)
		Inc(f[i],f[i-1]);
} 

int inv[K],ifc[K];

void init_fac(){
	inv[1]=ifc[1]=ifc[0]=1;
	for(int re i=2;i<K;++i){
		inv[i]=mul(mod-mod/i,inv[mod%i]);
		ifc[i]=mul(ifc[i-1],inv[i]);
	}
}
inline int get_f(ll x){return x<N?f[x]:F[n/x];}
void init_F(){
	for(int re i=n/N;i;--i){
		ll n=::n/i;F[i]=S3(n);
		for(ll re l=2,r,u;(u=n/l);l=r+1)
			Dec(F[i],mul(get_f(u),dec(S2(r=n/u),S2(l-1))));
	}
}

int y[K];
void init_lagrange(){
	y[0]=mod-1;k=tim+2;
	for(int re i=1;i<=k;++i)
		y[i]=add(y[i-1],po(i,tim));
}

int pr[K],sf[K]; 
int lagrange(ll n){
	int x=MD(n),res=0;pr[0]=sf[k+1]=1;
	for(int re i=1;i<=k;++i)
		pr[i]=mul(pr[i-1],dec(x,i));
	for(int re i=k;i;--i)
		sf[i]=mul(sf[i+1],dec(x,i));
	for(int re i=1;i<=k;++i){
		int coef=mul(ifc[i-1],ifc[k-i]);
		if((k^i)&1)coef=mod-coef;
		Inc(res,mul(mul(y[i],coef),mul(pr[i-1],sf[i+1])));
	}return res;
}

int g1[N],g2[N];
void min_25_sieve(){
	nn=sqrt(n);g1[0]=mod-1;
	for(int re i=1;i<=nn;++i){
		g1[i]=add(g1[i-1],po(i,tim));
		g2[i]=lagrange(n/i);
	}for(int re i=1;p[i]<=nn;++i){
		int p=::p[i],_g=g1[p-1],pw=po(p,tim);
		int w1=nn/p,w2=std::min((ll)nn,n/p/p);
		for(int re i=1;i<=w1;++i)
			Dec(g2[i],mul(pw,dec(g2[i*p],_g)));
		for(int re i=w1+1;i<=w2;++i)
			Dec(g2[i],mul(pw,dec(g1[n/i/p],_g)));
		for(int re i=nn;i>=(ll)p*p;--i)
			Dec(g1[i],mul(pw,dec(g1[i/p],_g)));
	}g1[0]=0;
}
inline int get_g(ll x){
	return x<=nn?g1[x]:g2[n/x];
}

void Main(){
	linear_sieves();init_fac();
	int T;scanf("%d",&T);
	while(T--){
		scanf("%lld%d",&n,&tim);
		init_F();++tim;
		init_lagrange();
		min_25_sieve();int ans=0;
		for(ll l=1,r,u;(u=n/l);l=r+1)
			Inc(ans,mul(dec(get_g(r=n/u),get_g(l-1)),get_f(u)));
		cout<<ans<<"\n";
	}
}

inline void file(){
#ifdef zxyoi
	freopen("math.in","r",stdin);
#endif
}
signed main(){file();Main();return 0;}