记\(f(x)=\)\(x\)的次大因数，那么\(sgcd(i,j)=f(gcd(i,j))\)。 下面来推式子：\[ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^nsgcd(i,j)^k&=\sum_{i=1}^n\sum_{j=1}^nf(gcd(i,j))^k&记d=gcd(i,j)\\ &=\sum_{d=1}^nf(d)^k\sum_{\frac i d=1}^{\lfloor \frac n d\rfloor}\sum_{\frac j d=1}^{\lfloor \frac n d\rfloor}[gcd(\frac i d,\frac j d)=1]&f(1)=0\\ &=\sum_{d=2}^nf(d)^k\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac n d\rfloor}[gcd(i,j)=1]\\ &=\sum_{d=2}^nf(d)^k(2\sum_{i=1}^{\lfloor \frac n d\rfloor}\varphi(i)-1) \end{aligned} \] 最后一步的括号是用欧拉函数的定义直接替换出来的。 我们发现，可以按\(\lfloor \frac n d \rfloor\)的取值数论分段，因为括号显然只受\(\lfloor \frac n d\rfloor\)的取值影响，关键是如何求\(f(x)^k\)的前缀和？ 记\(S_x=\sum_{i=2}^xf(x)^k\)， 令min_25中的\(g\)数组以\(s(x)=x^k\)的定义计算\(g\)，考虑由\(g_{i,j-1}\)推到\(g_{i,j}\)的时候，减去的是什么？是\(s(p_j)*(g_{\lfloor\frac n {p_j}\rfloor,j-1}-g_{p_j-1,j-1})\)，后面括号的部分是什么呢？恰好是那些最小质因子为\(p_j\)且除去\(p_j\)后剩余部分最小质因数不小于\(p_j\)的数的\(k\)次方之和，我们发现这些数的\(f\)值之和就是后面的括号。又因为每一个数至多被如此枚举到1次，所以对于每一个\(i\)，我们把括号的值累加起来，这就是那些\([2,i]\)中非质数的\(f\)值之和；再加上\([2,i]\)中的质数的\(f\)值之和，也就是质数个数，我们就求得了\(S_i\)。 min_25真牛逼。 所以就做完了，\(g\)的初始化需要用到自然数幂求和，考虑到\(k\)比较小，可以用第二类斯特林数求解。\[ \begin{aligned} Sum_k(n)&=\sum_{i=1}^ni^k\\ &=\sum_{j=1}^k\begin{Bmatrix}k\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1} \end{aligned} \]
#include <cstdio> #include <cmath> #include <algorithm> #include <map> using namespace std; typedef long long ll; typedef unsigned int ui; typedef map<ll,ui> mlu; const int SQRTN=32000,LEN=1000001; bool vis[LEN]; ui p[LEN],pcnt,pphi[LEN],s[51][51],pk[LEN]; ll n,k,sqrtn,m; ll a[SQRTN*2],cnt,pos1[SQRTN],pos2[SQRTN]; ui g[SQRTN*2],sum[SQRTN*2],g1[SQRTN*2]; mlu record; ui ksm(ui x,int y){ ui res=1; for(;y;x=x*x,y>>=1) if(y&1) res=res*x; return res; } void prework(){ pphi[1]=1; for(int i=2;i<LEN;i++){ if(!vis[i]){ p[++pcnt]=i; pphi[i]=i-1; } for(int j=1;j<=pcnt&&i*p[j]<LEN;j++){ int x=i*p[j]; vis[x]=true; if(i%p[j]==0){ pphi[x]=pphi[i]*p[j]; break; } pphi[x]=pphi[i]*(p[j]-1); } } for(int i=2;i<LEN;i++) pphi[i]+=pphi[i-1]; s[0][0]=s[0][1]=1; for(int i=1;i<=50;i++){ s[i][1]=1; for(int j=2;j<=i;j++) s[i][j]=s[i-1][j]*(ui)j+s[i-1][j-1]; } } inline int gp(ll x){ return x<=sqrtn?pos1[x]:pos2[n/x]; } void Diz(){ for(ll i=1,j;i<=n;i=j+1){ j=n/(n/i); a[++cnt]=n/i; } reverse(a+1,a+1+cnt); for(int i=1;i<=cnt;i++) if(a[i]<=sqrtn) pos1[a[i]]=i; else pos2[n/a[i]]=i; } ui sumup(ll l,ll r){ ll a=l+r,b=r-l+1; if(a&1) b/=2; else a/=2; return (ui)a*(ui)b; } ui getPhi(ll n){ if(n<LEN) return pphi[n]; mlu::iterator pos=record.find(n); if(pos!=record.end()) return pos->second; ui res=sumup(1,n); for(ll i=2,j;i<=n;i=j+1){ j=n/(n/i); res-=getPhi(n/i)*(j-i+1); } record[n]=res; return res; } ui getSumk(ll n){ ui res=0,t; for(ll j=1;j<=k;j++){ t=1; int md=(n-j+1)%(j+1); for(ll i=n-j+1;i<=n+1;i++,md++) if(!md||md==(j+1)) t*=(ui)i/(j+1); else t*=(ui)i; res+=t*s[k][j]; } return res; } void calc_g(){ for(int i=2;i<=cnt;i++) g[i]=getSumk(a[i])-1,g1[i]=a[i]-1,sum[i]=0; for(int j=1;j<=m;j++) for(int i=cnt;i>=2&&a[i]>=p[j]*p[j];i--){ g1[i]-=g1[gp(a[i]/p[j])]-g1[gp(p[j]-1)]; ui delta=(g[gp(a[i]/p[j])]-g[gp(p[j]-1)]); g[i]-=pk[j]*delta; sum[i]+=delta; } } ui calc(ll n){ return sum[gp(n)]+g1[gp(n)]; } int main prework(); scanf("%lld%lld",&n,&k); sqrtn=(ll)sqrt(n); m=upper_bound(p+1,p+1+pcnt,sqrtn)-p-1; for(int i=1;i<=m;i++) pk[i]=ksm(p[i],k); Diz(); calc_g(); ui ans=0,last=0,tmp; for(ll i=2,j;i<=n;i=j+1){ j=n/(n/i); tmp=calc(j); ans+=(getPhi(n/i)*2-1)*(tmp-last); last=tmp; } printf("%lu\n",ans); return 0; }
转载于:https://www.cnblogs.com/RogerDTZ/p/9196482.html
相关资源：数据结构—成绩单生成器