51nod 1847 奇怪的数学题

求解${\sum }_{i=1}^n{\sum }_{j=1}^{n}sgcd(i,j),sgcd$表示次大公约数,$n\le 10^{10}$

那么首先有$sgcd(i,j)=gcd(i,j)/mn(gcd(i,j))$就是gcd除以最小质因子。现在可以考虑交换枚举顺序先枚举gcd，然后最终得到的式子是这个形式${\sum }_{d=1}^n\frac{d^k}{mn(d)}({\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }2\phi (i)-1)$

然后看到这个形式我们就可以整除分块了，后面这一部分可以杜教筛，关键在于前面，这个东西和最小质因子相关，我们可以利用min25筛的第一部分，因为我们每次转移就相当于是枚举了最小质因子，那么剩下的就是答案，然后初始化的时候需要利用第二类斯特林数求自然数幂之和，因为模数不是质数。

#include<bits/stdc++.h>
#define LL long long
#define uint unsigned int
#define V inline void
#define I inline uint
#define FOR(i,a,b) for(register int i=a,end##i=b;i<=end##i;++i)
#define REP(i,a,b) for(register int i=a,end##i=b;i>=end##i;--i)
#define go(i,x) for(int i=hed[x];i;i=e[i].pre)
using namespace std;
inline int read()
{char x='\0';int fh=1,sum=0;for(x=getchar();x<'0'||x>'9';x=getchar())if(x=='-')fh=-1;for(;x>='0'&&x<='9';x=getchar())sum=sum*10+x-'0';return fh*sum;
}
const int N=1e6+9;
uint n,k;
I ksm(uint a,uint b)
{uint sum=1;while(b){if(b&1)sum*=a;b>>=1;a*=a;}return sum;
}
uint s[59][59];
I getsum(uint x)
{uint ret=0;FOR(i,1,k){uint sum=1;FOR(j,0,i){if((x+1-j)%(i+1)==0)sum*=(x+1-j)/(i+1);else sum*=x+1-j;}ret+=s[k][i]*sum;}return ret;
}
uint phi[N],pri[N],pcnt;
bool vis[N];
uint sp0[N],spk[N];
V oula(uint mx)
{vis[0]=vis[1]=true;phi[1]=1;for(uint i=2;i<=mx;i++){if(!vis[i])pri[++pcnt]=i,phi[i]=i-1,sp0[pcnt]=sp0[pcnt-1]+1,spk[pcnt]=spk[pcnt-1]+ksm(i,k);for(uint j=1;j<=pcnt&&i*pri[j]<=mx;j++){vis[i*pri[j]]=true;phi[i*pri[j]]=phi[i]*(pri[j]-1);if(i%pri[j]==0){phi[i*pri[j]]=pri[j]*phi[i];break;}}}FOR(i,1,mx)phi[i]+=phi[i-1];
}
map<uint,uint>sphi;
I getphi(uint x)
{if(x<=N-9)return phi[x];if(sphi.find(x)!=sphi.end())return sphi[x];uint ans=x*(x+1)/2;for(uint lp=2,rp=0;lp<=x;lp=rp+1){uint now=x/lp;rp=x/now;ans-=(rp-lp+1)*getphi(now);}sphi[x]=ans;return ans;
}
uint sq;
uint w[N],tot,ind1[N],ind2[N];
uint g0[N],gk[N],ans[N];
int main()
{n=read(),k=read();sq=ceil(sqrt(n));oula(1e6);s[0][0]=1;FOR(i,1,k)FOR(j,1,i)s[i][j]=j*s[i-1][j]+s[i-1][j-1];for(uint lp=1,rp=0;lp<=n;lp=rp+1){uint now=n/lp;rp=n/now;w[++tot]=now;g0[tot]=now-1,gk[tot]=getsum(now)-1;
//		cout<<now<<' '<<g0[tot]<<' '<<gk[tot]<<endl;//if(now<=sq)ind1[now]=tot;else ind2[n/now]=tot;}for(uint i=1;i<=pcnt&&pri[i]<=sq;i++){for(uint j=1;j<=tot&&1LL*pri[i]*pri[i]<=w[j];j++){uint now=w[j]/pri[i];uint h=(now<=sq)?ind1[now]:ind2[n/now];gk[j]-=ksm(pri[i],k)*(gk[h]-spk[i-1]);ans[j]+=gk[h]-spk[i-1];g0[j]-=g0[h]-sp0[i-1];
//			cout<<pri[i]<<" "<<w[j]<<' '<<g0[j]<<' '<<gk[j]<<' '<<ans[j]<<endl;}}uint sum=0;for(uint lp=1,rp=0;lp<=n;lp=rp+1){uint now=n/lp;rp=n/now;uint x=(rp<=sq)?ind1[rp]:ind2[n/rp];uint y=(lp-1<=sq)?ind1[lp-1]:ind2[n/(lp-1)];
//		cout<<"now "<<now<<' '<<getphi(now)<<endl;
//		cout<<"ans "<<ans[x]<<' '<<ans[y]<<endl;//
//		cout<<"g "<<g0[x]<<' '<<g0[y]<<endl;//sum+=(2*getphi(now)-1)*(ans[x]+g0[x]-ans[y]-g0[y]);
//		cout<<"sum "<<sum<<endl;//}printf("%u",sum);return 0;
}