【51nod1847】奇怪的数学题(Min_25筛+杜教筛)
题面
题解
这题有毒……不知为啥的错误调了半天……
令\(f(i)={sgcd(i)}\),那么容易看出\(f(i)\)就是\(i\)的次大质因子,用\(i\)除以它的最小质因子即可计算
于是题目所求即为
\[\sum_{i=1}^n\sum_{j=1}^n{f(\gcd(i,j))}^k
\]
\[\sum_{d=1}^n {f(d)}^k\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]
\]
\[\sum_{d=1}^n {f(d)}^k(2\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\varphi(i)-1)
\]
后面那个根据\(\varphi\)的定义就可以推出来
因为后面括号中的式子只与\({\left\lfloor\frac{n}{d}\right\rfloor}\)的值有关,我们可以考虑整除分块,那么括号里的东西就是一个杜教筛
那么我们现在的问题就是对于前半部分怎么快速求和了。我们考虑\(Min\_25\)筛的过程,其中\(g(n,j)\)表示所有\(1\)到\(n\)中的质数或最小质因子大于\(P_j\)的所有数的\(k\)次方之和,即
\[g(n,j)=\sum_{i=1}^ni^k[i\in P\ or\ min_p(i)>P_j]
\]
中间的转移中有这么一步
\[g(n,i)=g(n,j-1)-{P_j}^k(g({\left\lfloor\frac{n}{P_j}\right\rfloor},j-1)-\sum_{i=1}^{j-1}{P_i}^k)
\]
考虑后面减去的东西,是所有最小质因子等于\(P_j\)且自身为合数的数的\(k\)次方之和,即
\[\sum_xx^k[x\notin P\ and\ min_p(x)=P_j]
\]
然后我们惊喜地发现,后面减去的东西中,括号里的东西就是上面所有满足条件的\(x\)的\({f(x)}^k\)之和!
那么我们就可以直接\(Min\_25\)筛来计算\({f(x)}^k\)的前缀和了,那么一段区间只要差分一下就可以了
最后有个尴尬的问题,就是这个模数不是质数,所以我们在初始化的时候需要快速计算\(\sum_{i=1}^n i^k\)就不能拉格朗日插值了
那么只好用第二类斯特林数优化了
\[\begin{aligned}
\sum\limits_{i=1}^{n}i^K&=\sum_{i=1}^n\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}i^\underline{j}\\
&=\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}\sum_{i=1}^ni^\underline{j}\\
&=\sum_{j=1}^K\begin{Bmatrix}K\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1}
\end{aligned}\]
然后没有然后了
//minamoto
#include<bits/stdc++.h>
#define R register
#define int unsigned int
#define IT map<int,int>::iterator
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=1e5+5;
int w[N<<1],g[N<<1],h[N<<1],G[N<<1],id1[N],id2[N],p[N],pw[N],sp[N],phi[N],S[55][55];
map<int,int>mp;bitset<N>vis;int n,m,k,tot,sqr,id,ans,now,las;
int ksm(R int x,R int y){
R int res=1;
for(;y;y>>=1,x*=x)if(y&1)res*=x;
return res;
}
void init(int n){
phi[1]=1;
fp(i,2,n){
if(!vis[i])p[++tot]=i,pw[tot]=ksm(i,k),sp[tot]=sp[tot-1]+pw[tot],phi[i]=i-1;
for(R int j=1;j<=tot&&1ll*i*p[j]<=n;++j){
vis[i*p[j]]=1;
if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
fp(i,1,n)phi[i]+=phi[i-1];
S[0][0]=1;
fp(i,1,k){
S[i][1]=1;
fp(j,2,i)S[i][j]=S[i-1][j]*j+S[i-1][j-1];
}
}
int Phi(int n){
if(n<=sqr)return phi[n];
IT it=mp.find(n);
if(it!=mp.end())return it->second;
int res=n*(n+1)>>1;
for(R int i=2,j;i<=n;i=j+1)
j=n/(n/i),res-=Phi(n/i)*(j-i+1);
return mp[n]=res;
}
int calc(int n){
int res=0,tmp;
fp(i,1,min(k,n)){
tmp=S[k][i];
fp(j,n-i+1,n+1)(j%(i+1)==0)?tmp*=j/(i+1):tmp*=j;
res+=tmp;
}return res;
}
signed main(){
// freopen("testdata.in","r",stdin);
scanf("%u%u",&n,&k),init(sqr=sqrt(n));
for(R int i=1,j;i<=n;i=j+1){
j=n/(n/i),w[++m]=n/i;
w[m]<=sqr?id1[w[m]]=m:id2[n/w[m]]=m;
g[m]=calc(w[m])-1,h[m]=w[m]-1;
}
fp(j,1,tot)for(R int i=1;1ll*p[j]*p[j]<=w[i];++i){
id=(w[i]/p[j]<=sqr)?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
g[i]-=pw[j]*(g[id]-sp[j-1]);
h[i]-=h[id]-j+1;
G[i]+=g[id]-sp[j-1];
}
for(R int i=2,j;i<=n;i=j+1){
j=n/(n/i),id=(j<=sqr)?id1[j]:id2[n/j];
now=G[id]+h[id],ans+=((Phi(n/i)<<1)-1)*(now-las),las=now;
}
printf("%u\n",ans);
return 0;
}
深深地明白自己的弱小