题解：P15398 前缀和与差分 / coprime

传送门

首先先推个式子，考虑枚举序列中所有数的最大公因数，于是可以得到以下式子：

\[f(N,i)=\sum_{d=1}^i d^2 \sum_{a_1=1}^{\lfloor\frac{i}{d}\rfloor} \sum_{a_2=1}^{\lfloor\frac{i}{d}\rfloor} \dots \sum_{a_N=1}^{\lfloor\frac{i}{d}\rfloor} [\gcd(a_1,a_2,\dots,a_N)=1] \]

考虑莫反，可以推出以下式子：

\[\begin{aligned} f(N,i)&=\sum_{d=1}^i d^2 \sum_{d'=1}^{\lfloor\frac{i}{d}\rfloor} \sum_{a_1=1}^{\lfloor\frac{i}{d\times d'}\rfloor} \sum_{a_2=1}^{\lfloor\frac{i}{d\times d'}\rfloor} \dots \sum_{a_N=1}^{\lfloor\frac{i}{d\times d'}\rfloor} \mu(d')\\ &=\sum_{d=1}^i d^2 \sum_{d'=1}^{\lfloor\frac{i}{d}\rfloor} \mu(d') ({\lfloor\frac{i}{d\times d'}\rfloor})^N \end{aligned}\]

考虑差分，对于 \(i\)，求出 \(i\) 相对于 \(i-1\) 多出的贡献，记其为 \(s_i\)，则有：

\[s_i=\sum_{d|i} \sum_{d'|\frac{i}{d}} d^2 \mu(d') (({\frac{i}{d\times d'}})^N-({\frac{i}{d\times d'}}-1)^N) \]

考虑分步计算，先计算出 \(g_i=\sum_{d|i} d^2\mu(\frac{i}{d})\)，然后再计算出 \(s_i=\sum_{d|i} g_i\times (({\frac{i}{d\times d'}})^N-({\frac{i}{d\times d'}}-1)^N)\)。

注意到 \(f(i)=i^2\) 和 \(f(i)=\mu(i)\) 都是积性函数，于是可以狄利克雷卷积快速计算式子。

至于计算 \((({\frac{i}{d\times d'}})^N-({\frac{i}{d\times d'}}-1)^N)\)，可以先对于所有 \(1\le x\le M\) 求出所有 \(x^N\)，然后作差分，计算 \(x^N\) 可以线性筛做到线性，具体地，在质数处暴力计算，然后枚举合数时直接合并。

CODE：

#include<bits/stdc++.h>
#define ll long long
#define uint unsigned int
#define ull unsigned long long
using namespace std;
int n,m;
int tot=0,a[2000010];
bool b[20000010];
uint p[20000010],s[20000010],mu[20000010],v[20000010];
uint ss[20000010],ans[20000010];
const ll mod=998244353,mod1=1e9+7;
inline ll add(ll x,ll y){if(x+y>=mod) return x+y-mod;return x+y;}
inline ll add1(ll x,ll y){if(x+y>=mod1) return x+y-mod1;return x+y;}
inline ll red(ll x,ll y){if(x-y<0) return x-y+mod;return x-y;}
inline ll qpow(ll x,ll y){ll sum=1;while(y){if(y&1) sum*=x,sum%=mod;x*=x,x%=mod,y>>=1;}return sum;}
inline void sss()
{
	b[0]=b[1]=1,mu[1]=1,v[1]=1,ss[1]=1;
	for(register int i=2;i<=m;i++)
	{
		if(!b[i]) tot++,a[tot]=i,v[i]=qpow(i,n),mu[i]=mod-1,p[i]=i,s[i]=1;
		for(register int j=1;j<=tot&&1ull*i*a[j]<=m;j++)
		{
			int x=i*a[j];
			b[x]=1,v[x]=1ull*v[i]*v[a[j]]%mod,p[x]=a[j];
			if(!(i%a[j])){s[x]=s[i];break;}
			else s[x]=i;
			mu[x]=1ull*(mod-1)*mu[i]%mod;
		}
		if(s[i]==1)
		{
			ull sssss=1;
			while(sssss<=i) ss[i]=add(ss[i],1ull*sssss*sssss%mod*mu[i/sssss]%mod),sssss*=p[i];
		}
		else ss[i]=1ull*ss[s[i]]*ss[i/s[i]]%mod;
	}
}
int main()
{
	scanf("%d%d",&n,&m);
	sss();
	for(register int i=m;i>=1;i--) v[i]=red(v[i],v[i-1]);
	for(register int i=1;i<=m;i++) ans[i]=v[i];
	for(register int i=1;i<=tot;i++)
		for(register int j=m/a[i];j>=1;j--)
		{
			ull sssss=a[i];int x=a[i]*j;
			while(!(x%sssss)) ans[x]=add(ans[x],1ull*ss[sssss]*ans[x/sssss]%mod),sssss*=a[i];
		}
	for(register int i=1;i<=m;i++) ans[i]=add(ans[i],ans[i-1]);
	ull S=0,ssum=1;
	for(register int i=m;i>=1;i--) S=add1(S,ssum*ans[i]%mod1),ssum=ssum*mod%mod1;
	printf("%lld\n",S);
	return 0;
}