【2019 Multi-University 3 E】Easy Math Problem 题解

题目大意

  求：

$$\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k lcm(i,j) [\gcd\ is\ prime]$$

  $1 \leq n \leq 10^{10},~1 \leq k \leq 100$
  多测，$T \leq 5$，时限 10s

$$\\$$ $$\\$$ $$\\$$

题解

$$\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k lcm(i,j) ~[\gcd~is~prime] \\ =\sum_{d=1}^n d^{k+1}[d~is~prime] \sum_{i'=1}^{\lfloor \frac nd \rfloor} \sum_{j'=1}^{\lfloor \frac nd \rfloor} i'j'[\gcd(i,j)=1] \\$$

  后面那玩意儿是个经典转换：

$$\begin{aligned} & \sum_{i=1}^m \sum_{j=1}^m ij[\gcd(i,j)=1] \\ =& \bigg( 2 \times \sum_{i=2}^m i \sum_{j=1}^{i-1} j[\gcd(i,j)=1] \bigg) + \bigg( \sum_{i=1}^m i^2[\gcd(i,i)=1] \bigg)\\ =& \bigg( 2 \times \sum_{i=2}^m i \frac{i \phi(i)}2 \bigg) +1 \\ =& \sum_{i=1}^m i^2 \phi(i) \end{aligned}$$

  因此原式等于

$$\sum_{d=1}^n d^{k+1} [d~is~prime] \sum_{i=1}^{\lfloor \frac nd \rfloor} i^2 \phi(i)$$

  于是对 $\lfloor \frac nd \rfloor$ 分块，前面求质数幂和用洲阁筛或 min25，后面求 $\sum i^2\phi(i)$ 用杜教筛或 min25。

扯淡

  symbol 推法（可参考 WC2016 第二课堂课件）居然不是万能的。。。

$$\begin{aligned} &\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k lcm(i,j) ~[\gcd~is~prime] \\ =&\sum_{d=1}^n d^{k+1}[d~is~prime] \sum_{i'=1}^{\lfloor \frac nd \rfloor} \sum_{j'=1}^{\lfloor \frac nd \rfloor} i'j'[\gcd(i,j)=1] \\ =&\sum_{d=1}^n d^{k+1}[d~is~prime] \sum_{D=1}^{\lfloor \frac nd \rfloor} \mu(D)D^2\frac{(1+\lfloor \frac n{dD} \rfloor)^2 (\lfloor \frac n{dD} \rfloor)^2}4 \\ =&\sum_{T=1}^n \frac{(1+\lfloor \frac nT \rfloor)^2 (\lfloor \frac nT \rfloor)^2}4 T \sum_{d|T,~d~is~prime} d^{k-1}\mu(\frac Td) \end{aligned}$$

  以前是觉得，这种变量前移的推法能通杀一切，即便跟标准做法不一样，但最后会推出来相同的东西，或者我这样推也是能做的（甚至更快）。做这题之前，这个想法都是对的。
  直到遇见了这题。。。后面那坨居然不能预处理又不能筛？？

代码

#include<bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;

typedef long long LL;

const int maxsqrtn=2e6+5, maxk=105;
const LL mo=1e9+7;

LL n;
int sqrtn,k;

LL mi(LL x,LL y)
{
	LL re=1;
	for(; y; y>>=1, x=x*x%mo) if (y&1) re=re*x%mo;
	return re;
}

LL phi[maxsqrtn],Sp[maxsqrtn],inv[maxsqrtn];
int p0,p[maxsqrtn],Np[maxsqrtn];
bool bz[maxsqrtn];
void Pre(int n)
{
	p0=0;
	memset(bz,0,sizeof(bz));
	memset(Np,0,sizeof(Np));
	phi[1]=1;
	fo(i,2,n)
	{
		if (!bz[i])
		{
			p[++p0]=i;
			phi[i]=i-1;
			Np[i]=1;
			Sp[p0]=(Sp[p0-1]+mi(i,k))%mo;
		}
		fo(j,1,p0)
		{
			if ((LL)i*p[j]>n) break;
			bz[i*p[j]]=1;
			if (i%p[j]==0)
			{
				phi[i*p[j]]=phi[i]*p[j];
				break;
			} else phi[i*p[j]]=phi[i]*(p[j]-1);
		}
	}
	inv[1]=1;
	fo(i,2,n)
	{
		phi[i]=(phi[i-1]+phi[i]*i%mo*i%mo)%mo;
		Np[i]+=Np[i-1];
		inv[i]=(-(mo/i)*inv[mo%i])%mo;
	}
}

LL x[maxk],y[maxk],w[maxk],invw[maxk];
void lagrange_add(int n)
{
	fo(i,1,n-1) (w[i]*=(x[i]-x[n]))%=mo;
	w[n]=1;
	fo(i,1,n-1) (w[n]*=(x[n]-x[i]))%=mo;
}
LL lagrange_get(int n,LL nx)
{
	if (nx<=n) return y[nx];
	nx%=mo;
	LL sum=0, l=1;
	fo(i,1,n)
	{
		(l*=(nx-x[i]))%=mo;
		(sum+=y[i]*invw[i]%mo*(nx-x[i]<=maxsqrtn-5 ?inv[nx-x[i]] :mi(nx-x[i],mo-2)))%=mo;
	}
	return l*sum%mo;
}

LL mw[maxsqrtn],g[maxsqrtn],id1[maxsqrtn],id2[maxsqrtn];
int w0;
LL min25_g(LL n)
{
	w0=0;
	for(LL i=1, j; i<=n; i=j+1)
	{
		j=n/(n/i);
		mw[++w0]=n/i;
		if (mw[w0]<=sqrtn) id1[mw[w0]]=w0; else id2[j]=w0;
		g[w0]=lagrange_get(k+2,mw[w0])-1;
	}
	fo(j,1,Np[sqrtn])
		for(int i=1; i<=w0 && (LL)p[j]*p[j]<=mw[i]; i++)
		{
			int id=(mw[i]/p[j]<=sqrtn) ?id1[mw[i]/p[j]] :id2[n/(mw[i]/p[j])];
			(g[i]-=(Sp[j]-Sp[j-1])*(g[id]-Sp[j-1]))%=mo;
		}
}

unordered_map<LL,LL> f;
LL inv2,inv6;
LL sum(LL n,int id)
{
	n%=mo;
	if (id==2) return n*(n+1)%mo*inv6%mo*(2*n+1)%mo;
		else if (id==3)
		{
			LL t=n*(n+1)%mo*inv2%mo;
			return t*t%mo;
		}
}
LL Sphi(LL x,int id)
{
	if (x<=maxsqrtn-5) return phi[x];
	if (f.count(x)) return f[x];
	
	LL re=sum(x,id+1);
	for(LL i=2, j; i<=x; i=j+1)
	{
		j=x/(x/i);
		(re-=(sum(j,id)-sum(i-1,id)+mo)%mo*Sphi(x/i,id)%mo)%=mo;
	}
	(re+=mo)%=mo;
	
	return f[x]=re;
}

int T;
int main()
{
	inv2=mi(2,mo-2);
	inv6=mi(6,mo-2);
	
	scanf("%d",&T);
	while (T--)
	{
		scanf("%lld %d",&n,&k); k++;
		sqrtn=sqrt(n);
		
		Pre(maxsqrtn-5);
		
		fo(i,1,k+2)
		{
			x[i]=i, y[i]=(y[i-1]+mi(i,k))%mo;
			lagrange_add(i);
		}
		fo(i,1,k+2) invw[i]=mi(w[i],mo-2);
		
		min25_g(n);
		
		LL ans=0;
		for(LL i=1, j; i<=n; i=j+1)
		{
			j=n/(n/i);
			int idi=((i-1)<=sqrtn) ?id1[(i-1)] :id2[n/(i-1)] ;
			int idj=(j<=sqrtn) ?id1[j] :id2[n/j] ;
			(ans+=(g[idj]-g[idi]+mo)%mo*Sphi(n/i,2))%=mo;
		}
		
		printf("%lld\n",(ans+mo)%mo);
	}
}