【51NOD1847】奇怪的数学题（杜教筛，min_25筛，第二类斯特林数解决自然数幂求和）

设 $f(n)$ 表示 $n$ 的次大因数。
$$\begin{array}{rl}& \sum _{i=1}^n\sum _{j=1}^{n}f(\gcd (i,j){)}^k\\ =& \sum _{d=1}^{n}f(d{)}^k\sum _{i=1}^{(n/d)}\sum _{j=1}^{(n/d)}[\gcd (i,j)=1]\\ =& \sum _{d=1}^{n}f(d{)}^k\left(2\sum _{i=1}^{(n/d)}\phi (i)-1\right)\end{array}$$
（其中第二步推到第三步是根据 $\phi$ 的定义，括号内 $-1$ 是因为 $i=j=1$ 的情况我们算了两次）

我们数论分块，然后后面括号里的 $\phi$ 的前缀和可以用杜教筛解决。

对于前面的 $f(d{)}^k$ 的前缀和，我们分为质数和合数两部分。

对于质数，$f(p)=1$，于是询问的就是区间质数个数，这个可以用 min_25 筛解决。

对于合数，观察 min_25 筛中推导 $g_k(n,i)=\sum _{x=1}^n[x\in \mathbb{P}\vee \mathrm{lpf}(x)>p_i]x^k$ 的过程：$g_k(n,i)$ 由 $g_k(n,i-1)$ 减去最小质因数为 $p_i$ 的合数的贡献得到。所以我们只需要加上求 $g_k(n,i)$ 过程中用 $g_k(n,i-1)$ 减去的东西即可。

但是注意初始化 $g_k(n,0)=\sum _{x=2}^n{x}^k$ 时，不能直接用拉格朗日插值来求自然数幂求和，因为模数 $2^{32}$ 不支持取逆。

那么就需要用到第二类斯特林数了：
$$\sum _{i=0}^n{i}^k=\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}\frac{(n+1{)}^{\underline{j+1}}}{j+1}$$
我们 $O(k^2)$ 预处理第二类斯特林数，然后每次 $O(k^2)$ 查询。

#include<bits/stdc++.h>

#define K 55
#define ll unsigned int

using namespace std;

inline ll poww(ll a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=ans*a;
		a=a*a;
		b>>=1;
	}
	return ans;
}

inline 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<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

int n,k;

namespace s1
{
	const int maxn=100000,N=maxn+10;
	const int SN=10010;
	int cnt,p[N];
	bool notprime[N];
	ll phi[N],sphi[SN];
	void init()
	{
		phi[1]=1;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		for(int i=2;i<=maxn;i++)
		{
			if(!notprime[i])
			{
				p[++cnt]=i;
				phi[i]=i-1;
			}
			for(int j=1,x;j<=cnt&&(x=i*p[j])<=maxn;j++)
			{
				notprime[x]=1;
				if(!(i%p[j]))
				{
					phi[x]=phi[i]*p[j];
					break;
				}
				phi[x]=phi[i]*phi[p[j]];
			}
		}
		for(int i=2;i<=maxn;i++)
			phi[i]=phi[i-1]+phi[i];
	}
	int get(int x)
	{
		return n/x;
	}
	ll S(int n)
	{
		if(n<=maxn) return phi[n];
		if(sphi[get(n)]) return sphi[get(n)];
		ll ans=1ll*n*(n+1)/2;
		for(int l=2,r;l<=n;l=r+1)
		{
			r=n/(n/l);
			ans-=S(n/l)*(r-l+1);
		}
		return sphi[get(n)]=ans;
	}
	ll query(int n)
	{
		return S(n);
	}
}

namespace s2
{
	const int N=32000;
	int sn,tot,b[N<<1],id1[N],id2[N];
	int cnt,p[N];
	bool notprime[N];
	ll s[K][K];
	ll powk[N],gpk[N];
	ll g0[N<<1],gk[N<<1];
	ll ans1[N<<1];
	void init()
	{
		for(int i=2;i<=sn;i++)
		{
			if(!notprime[i])
			{
				p[++cnt]=i,powk[cnt]=poww(i,k);
				gpk[cnt]=gpk[cnt-1]+powk[cnt];
			}
			for(int j=1,x;j<=cnt&&(x=i*p[j])<=sn;j++)
			{
				notprime[x]=1;
				if(!(i%p[j])) break;
			}
		}
	}
	ll getsum(int n)
	{
		ll ans=0;
		for(int j=0;j<=k;j++)
		{
			ll tmp=1;
			for(int l=0;l<=j;l++)
			{
				if((n+1-l)%(j+1)) tmp*=(n+1-l);
				else tmp*=(n+1-l)/(j+1);
			}
			ans+=s[k][j]*tmp;
		}
		return ans;
	}
	int get(int x)
	{
		return x<=sn?id1[x]:id2[n/x];
	}
	void work()
	{
		sn=sqrt(n);
		init();
		s[0][0]=1;
		for(int i=1;i<=k;i++)
			for(int j=1;j<=k;j++)
				s[i][j]=s[i-1][j-1]+s[i-1][j]*j;
		for(int l=1,r;l<=n;l=r+1)
		{
			r=n/(n/l);
			int x=n/l;
			b[++tot]=x;
			if(x<=sn) id1[x]=tot;
			else id2[n/x]=tot;
			g0[tot]=x-1;
			gk[tot]=getsum(x)-1;
		}
		for(int i=1;i<=cnt;i++)
		{
			int pi2=p[i]*p[i];
			for(int j=1;j<=tot&&pi2<=b[j];j++)
			{
				int v=get(b[j]/p[i]);
				g0[j]=g0[j]-(g0[v]-(i-1));
				gk[j]=gk[j]-powk[i]*(gk[v]-gpk[i-1]);
				ans1[j]+=gk[v]-gpk[i-1];
			}
		}
	}
	ll query(int n)
	{
		return ans1[get(n)]+g0[get(n)];
	}
}

int main()
{
	n=read(),k=read();
	s1::init();
	s2::work();
	ll ans=0,lst=0;
	for(int l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ll now=s2::query(r);
		ans+=(now-lst)*(s1::query(n/l)*2-1);
		lst=now;
	}
	cout<<ans<<endl;
	return 0;
}
/*
10 1
*/
/*
15 1
*/
/*
1000000 50
*/