[LOJ]#572. 「LibreOJ Round #11」Misaka Network 与求和 min_25筛+杜教筛

Solution

推一下式子，容易得到一个线性做法：$$\sum _{d=1}^n{f}^k(d)((2\sum _{i=1}^{\lfloor \frac{n}{i}\rfloor }\phi (i))-1)$$
这个东西数论分块加速一下，只需要快速求欧拉函数的前缀和和$f$的$k$次方的前缀和就可以解决了。欧拉函数前缀和直接杜教筛就搞定了，而$f$的$k$次方前缀和则需要用min_25筛。设$S(n,j)$表示$\le n$且最小质因子$\ge p_j$(第$j$个质数)的所有数$i$的$f^k(i)$之和。先默认$f(质数)=0$，最后再加上区间质数个数，考虑枚举最小的质因数及其次数进行转移，那么有：
当$p_k$不是次大质因数时，有$$S(n,j)+=\sum _{k\ge j}\sum _{e=1}^{{p_k}^{e+1}\le n}S(\lfloor \frac{n}{{p_k}^e}\rfloor ,k+1)$$
当$p_k$是次大质因数时，有$$S(n,j)+=\sum _{i=p_k}^{\lfloor \frac{n}{{p_k}^e}\rfloor }[i是质数]$$
区间内的质数个数也是min_25筛的基本应用。
由于有多次询问，需要对$S$的值进行记忆化，但其它信息不用重新处理，因为每次询问的值都是$\lfloor \frac{n}{i}\rfloor$的一个值。
复杂度我也不知道啊……

Code

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define pa pair<int,int>
#define ui unsigned int
const int Maxn=2000010;
const int inf=2147483647;
int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return x*f;
}
ui Pow(ui x,int y)
{
	if(!y)return 1;
	ui t=Pow(x,y>>1),re=t*t;
	if(y&1)re=re*x;
	return re;
}
int l,r,N,N1,n,k;
int tot=0,id1[62100],id2[62100];ui prime[Maxn],pw[Maxn],w[Maxn],g[Maxn];int lw=0;
ui phi[Maxn],sum[Maxn];
bool mark[Maxn];
int M;
void pre(int n)
{
	N1=(int)pow(n,2.0/3.0);
	N=(int)sqrt(n);phi[1]=1;int tot1;
	for(int i=2;i<=N1;i++)
	{
		if(!mark[i])
		{
			prime[++tot]=i,phi[i]=i-1;
			if(i<=N)pw[i]=Pow(i,k);
		}
		if(i==N)tot1=tot;
		for(int j=1;j<=tot&&prime[j]*i<=N1;j++)
		{
			mark[prime[j]*i]=true;
			if(i%prime[j]==0)
			{
				phi[prime[j]*i]=prime[j]*phi[i];
				break;
			}
			phi[prime[j]*i]=(prime[j]-1)*phi[i];
		}
	}
	sum[0]=0;for(int i=1;i<=N1;i++)sum[i]=sum[i-1]+phi[i];
	tot=tot1;
	int pos;lw=0;M=n;
	for(int i=1;i<=n;i=pos+1)
	{
		pos=n/(n/i);
		w[++lw]=n/i;
		if(n/i<=N)id1[n/i]=lw;
		else id2[n/(n/i)]=lw;
		g[lw]=w[lw]-1;
	}
	for(int j=1;j<=tot;j++)
	for(int i=1;i<=lw&&prime[j]*prime[j]<=w[i];i++)
	{
		int id;
		if(w[i]/prime[j]<=N)id=id1[w[i]/prime[j]];
		else id=id2[n/(w[i]/prime[j])];
		g[i]=g[i]-g[id]+j-1;
	}
} 
unordered_map<int,ui>s[62100];
ui S(int n,int m)
{
	if(n<2||prime[m]>n)return 0;
	if(s[m][n])return s[m][n];
	ui re=0;
	for(int j=m;j<=tot&&(LL)prime[j]*prime[j]<=n;j++)
	{
		int t=prime[j];
		for(int k=1;(LL)t*prime[j]<=n;k++)
		{
			re+=S(n/t,j+1);
			int id;
			if(n/t<=N)id=id1[n/t];
			else id=id2[M/(n/t)];
			re+=pw[prime[j]]*(g[id]-j+1);
			t*=prime[j];
		}
	}
	return s[m][n]=re;
}
int id3[62100],id4[62100],tot2=0;ui v[62100];
ui sp(int n)
{
	if(n<=N1)return sum[n];
	int id;
	if(n<=N)id=id3[n];else id=id4[M/n];
	if(v[id])return v[id];
	if(n<=N)id=id3[n]=++tot2;else id=id4[M/n]=++tot2;
	ui re;if(n&1)re=(n+1)/2*n;else re=n/2*(n+1);
	int pos;
	for(int i=2;i<=n;i=pos+1)
	{
		pos=n/(n/i);
		re-=sp(n/i)*(pos-i+1);
	}
	return v[id]=re;
}
int main()
{
	n=read(),k=read();
	pre(n);
	int pos;ui ans=0,lst=0;
	for(int i=1;i<=n;i=pos+1)
	{
		pos=n/(n/i);
		int id;
		if(pos<=N)id=id1[pos];else id=id2[n/pos];
		ui now=S(pos,1)+g[id];
		ans+=(now-lst)*(sp(n/i)*2-1);
		lst=now;
	}
	printf("%u",ans);
}