放松型休闲无脑小清新签到题

题目大意：

求\(\sum_{i=1}^n\sum_{j=1}^n [gcd(i,j)是质数]*gcd(i,j)^k,1\leq n\leq 10^9,1\leq k\leq 100\)

题解：

令\(f[x]=[x是质数]*x^k,F[x]=\sum_{i=1}^{x}f[i]\)
\(\sum_{i=1}^{n}\sum_{j=1}^{n}f[gcd(i,j)]=\sum_{d=1}^{n}f[d]\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=d]\)
\(=\sum_{d=1}^{n}f[d]\sum_{i=1}^{ \left\lfloor{\frac{n}{d}}\right\rfloor}\sum_{j=1}^{ \left\lfloor{\frac{n}{d}}\right\rfloor}e[gcd(i,j)]=\sum_{d=1}^{n}f[d]\sum_{i=1}^{ \left\lfloor{\frac{n}{d}}\right\rfloor}\sum_{j=1}^{ \left\lfloor{\frac{n}{d}}\right\rfloor}\sum_{e|i,e|j}^{}μ[e]\)
\(=\sum_{d=1}^{n}f[d]\sum_{e=1}^{ \left\lfloor{\frac{n}{d}}\right\rfloor}{ \left\lfloor{\frac{n}{de}}\right\rfloor}^{2}=\sum_{(t=de)=1}^{n}{ \left\lfloor{\frac{n}{t}}\right\rfloor}^{2}\sum_{de=t}^{}f[d]*μ[e]\)
分块，需要求\(G[x]=\sum_{t=1}^{x}\sum_{de=t}^{}f[d]*μ[e]=\sum_{e=1}^{x}μ[e]\sum_{d=1}^{ \left\lfloor{\frac{x}{e}}\right\rfloor}f[d]\sum_{e=1}^{x}μ[e]F[\left\lfloor{\frac{x}{e}}\right\rfloor]\)
再分块，需要求\(F[x]\)
发现和质数有关，考虑\(Min25\)筛。初始时需要自然数幂和，要伯努利数。
另：为卡常，可预处理部分\(F[]\)或不用\(Map\)记忆化杜教筛。

代码：

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
const int MOD=998244353;
int add(int a,int b)
{a+=b;return a>=MOD?a-MOD:a;}
int sub(int a,int b)
{a-=b;return a<0?a+MOD:a;}
int mul(int a,int b)
{return 1LL*a*b%MOD;}
int gcd(int a,int b)
{return b?gcd(b,a%b):a;}
int ksm(int a,int b)
{
	int ans=1;
	for(;b;b>>=1,a=mul(a,a))
		if(b&1)ans=mul(ans,a);
	return ans;
}
int n,k;
int s;
const int Q=1<<17;
int id1[Q],id2[Q],p[1<<18];
int ch[Q],mu[1<<20];
int tot=0;
int c[1<<7][1<<7];
int b[1<<7];
int f[Q];
int Gid(int x)
{return x<=s?id1[x]:id2[n/x];}
namespace Min25{
	int val[Q];
	int smf[Q];
	int inv[1<<7];
	int ch[1<<20]={0};
	void Prime(){
		int maxn=1e6;
		mu[1]=1;
		for(int i=2;i<=maxn;i++){
			if(!ch[i])p[++tot]=i,mu[i]=MOD-1;
			for(int j=1,t;j<=tot&&(t=i*p[j])<=maxn;j++){
				ch[t]=1;
				if(i%p[j]==0){
					mu[t]=0;
					break;
				}
				mu[t]=sub(0,mu[i]);
			}
		}
		for(int i=2;i<=maxn;i++)mu[i]=add(mu[i],mu[i-1]);
		for(int i=1;i<=tot;i++)
			if(p[i]>s){
				tot=i-1;
				return;
			}
	}
	void Bo(){
		for(int i=0;i<=123;i++){
			c[i][0]=1;
			for(int j=1;j<=i;j++)
				c[i][j]=add(c[i-1][j-1],c[i-1][j]);
		}
		b[0]=1;
		inv[1]=1;
		for(int i=1;i<=111;i++){
			int tmp=0;
			for(int j=0;j<i;j++)
				tmp=add(tmp,mul(c[i+1][j],b[j]));
			inv[i+1]=ksm(i+1,MOD-2);
			b[i]=mul(inv[i+1],sub(0,tmp));
		}
	}
	int Psum(int x){
		if(x<0){
			int als=0;
			for(int i=1;i<=x;i++)
				als=add(als,ksm(i,k));
			return als;
		}
		int als=0;
		x=(x+1)%MOD;
		for(int i=1,ow=x;i<=k+1;i++)
			als=add(als,mul(mul(c[k+1][i],b[k+1-i]),ow)),ow=mul(ow,x);
		return mul(inv[k+1],als);
	}
	void Init(){
		int maxn=0;
		Prime();Bo();
		for(int i=1,p1,p2;i<=n;i=p2+1)
		{
			p1=n/i,p2=n/p1;
			if(p1<=s)id1[p1]=++maxn;
			else id2[p2]=++maxn;
			val[maxn]=p1;
			f[maxn]=sub(Psum(p1),1);
		}
		for(int j=1;j<=tot;j++){
			int pr=p[j],owo=ksm(pr,k);
			smf[j]=add(smf[j-1],owo);
			for(int i=1;i<=maxn&&pr*pr<=val[i];i++)
				f[i]=sub(f[i],mul(owo,sub(f[Gid(val[i]/pr)],smf[j-1])));
		}
	}
}
int Res[1<<12];
int Gmu(int x)
{
	if(x<=1000000)return mu[x];
	if(Res[n/x])return Res[n/x]; 
	int als=0;
	for(int i=2,las;i<=x;i=las+1){
		las=x/(x/i);
		als=add(als,mul((las-i+1),Gmu(x/i)));
	}
	als=sub(1,als);
	return (Res[n/x]=als);
}
void Test(){
	int ppp;
	cin>>ppp;
	for(int i=1;i*i<=ppp;i++)
		if(ppp/(ppp/i)!=i)exit(0);
	Test();
}
int Gsum(int x)
{
	int als=0;
	for(int i=1,las;i<=x;i=las+1){
		las=x/(x/i);
		als=add(als,mul(sub(Gmu(las),Gmu(i-1)),f[Gid(x/i)]));
	}
	return als;
}
int main()
{
	cin>>n>>k;
	s=sqrt(n)+1;
	Min25::Init();
	int als=0;
	for(int i=1,las,lv=0,nv;i<=n;i=las+1){
		las=n/(n/i);
		nv=Gsum(las);
		als=add(als,mul(mul(n/i,n/i),sub(nv,lv)));
		lv=nv;
	}
	printf("%d",als);
	return 0;
}