Description
给出 $N,K$ ,请计算下面这个式子:
$sum _{i=1}^n sum_{j=1}^n sgcd(i,j)^k$
其中,$sgcd(i, j)$表示$(i, j)$的所有公约数中第二大的,特殊地,如果$gcd(i, j) = 1$, 那么$sgcd(i, j) = 0$。
考虑到答案太大,请输出答案对$2^{32}$取模的结果.
$1 leq N leq 10^9,1 leq K leq 50$
Solution
设$f(x)$为$x$的次大的因数，那么原式可化为
egin{align}
& sum _{i=1}^n sum_{j=1}^n sgcd(i,j)^k\
= & sum _{i=1}^n sum _{j=1}^n f((i,j))^k\
= & sum _{d=1}^n f^k(d) sum _{i=1}^{lfloor frac{n}{d}
floor }sum _{j=1}^{lfloor frac{n}{d}
floor}[(i,j)=1]\
= & sum _{d=1}^n f^k(d) left ( 2sum _{i=1}^{lfloor frac {n}{d}
floor }varphi(i)-1
ight )
end{align}
再定义$G_n=sum _{i=1}^n f(i)$，定义Min_25中的$g$函数为符合要求的$x^k$之和
发现在做$g$转移时：
$$g_{i,j}=g_{i-1,j}-p_j^k(g_{lfloor frac{i}{p_j} floor ,j-1 }-g_{p_{j-1},j-1})$$
括号内的内容即为该数本身除以它的最小质因数，就是题中的次大因数，所以在转移$g$的同时记录括号中的数的和，这样可以得到所有整除点的$G$值
$g$函数初始化的时候要求自然数幂之和，使用第二类斯特林数和普通幂的转化
欧拉函数可以用杜教筛得到整除点的值
#include<unordered_map> #include<iostream> #include<cstdio> #include<cmath> using namespace std; int tot,cnt,id1[100005],id2[100005]; unsigned int phi[1000005],S[60][60],pw[1000005],w[200005],g1[200005],g0[200005],sp[1000005],G[1000005],ans,pre,now; long long n,K,prime[1000005],sq; bool vst[1000005]; unordered_map<int,unsigned int>mp; inline int read(){ int f=1,w=0; char ch=0; while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();} while(ch>='0'&&ch<='9')w=(w<<1)+(w<<3)+ch-'0',ch=getchar(); return f*w; } unsigned int ksm(unsigned int a,unsigned int p){ unsigned int ret=1; while(p){ if(p&1)ret*=a; a*=a,p>>=1; } return ret; } unsigned int calc(long long x){ unsigned int ret=0,temp; for(unsigned int i=1;i<=K;i++){ int t=(x-i+1)%(i+1); temp=1; for(long long j=x-i+1;j<=x+1;j++,++t){ if(t>=i+1)t-=i+1; if(!t)temp*=(unsigned int)j/(i+1); else temp*=(unsigned int)j; } ret+=S[K][i]*temp; } return ret; } unsigned int djs(long long x){ if(x<=1000000)return phi[x]; if(mp.find(x)!=mp.end())return mp[x]; unsigned int ret=0; if(x&1)ret=(x+1)/2*x; else ret=x/2*(x+1); for(unsigned int i=2;i<=x;){ unsigned int j=x/(x/i); ret-=djs(x/i)*(j-i+1),i=j+1; } return mp[x]=ret; } int main(){ n=read(),K=read(),sq=sqrt(n),phi[1]=1; for(int i=2;i<=1000000;i++){ if(!vst[i])prime[++tot]=i,phi[i]=i-1; for(int j=1;j<=tot&&i*prime[j]<=1000000;j++){ vst[i*prime[j]]=true; if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j];break;} else phi[i*prime[j]]=phi[i]*phi[prime[j]]; } } for(int i=1;i<=1000000;i++)phi[i]+=phi[i-1]; for(int i=1;i<=K;i++){ S[i][1]=1; for(unsigned int j=2;j<=i;j++)S[i][j]=S[i-1][j]*j+S[i-1][j-1]; } for(int i=1;i<=tot;i++)pw[i]=ksm(prime[i],K),sp[i]=sp[i-1]+pw[i]; for(int i=1;i<=n;){ int j=n/(n/i); w[++cnt]=n/i,g1[cnt]=calc(w[cnt])-1,g0[cnt]=w[cnt]-1,n/i<=sq?id1[n/i]=cnt:id2[n/(n/i)]=cnt,i=j+1; } for(int i=1;i<=tot;i++)for(int j=1;j<=cnt&&1ll*prime[i]*prime[i]<=w[j];j++){ int id=w[j]/prime[i]<=sq?id1[w[j]/prime[i]]:id2[n/(w[j]/prime[i])]; g1[j]-=pw[i]*(g1[id]-sp[i-1]),g0[j]-=g0[id]-(i-1),G[j]+=g1[id]-sp[i-1]; // cout<<g0[j]<<" "; } for(int i=2;i<=n;){ int j=n/(n/i),id=j<=sq?id1[j]:id2[n/j]; now=G[id]+g0[id],ans+=(djs(n/i)*2-1)*(now-pre),pre=now,i=j+1; } printf("%u",ans); return 0; }//3656808373