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}^{?\frac{n}{d}?}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; }

總結

以上是生活随笔為你收集整理的51nod 1847 奇怪的数学题（数论/min25筛/杜教筛/斯特林数）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。