【Project Euler】530 GCD of Divisors 莫比烏斯反演
【題目】GCD of Divisors
【題意】給定f(n)=Σd|n gcd(d,n/d)的前綴和F(n),n=10^15。
【算法】莫比烏斯反演
【題解】參考:任之洲數論函數.pdf
這個範圍顯然杜教篩也是做不了的,而且考慮直接化簡f(n)也遇到了困難,所以考慮將前綴和的Σ一起化簡。
$$F(n)=\sum_{i=1}^{n}\sum_{d|i}(d,\frac{i}{d})$$
這一步很常見的是第一重改為枚舉倍數,但這樣化簡後面就推不下去了。
這道題必須最後轉成$\sigma_0(n)$才能解出來。
所以直接枚舉gcd值
$$F(n)=\sum_{d=1}^{n}d\sum_{i=1}^{n}\sum_{g|i}[(g,\frac{i}{g})=d]$$
這裏gcd(g,i/g)=d,說明i中必須至少包含2個d,那麽令i=i/d^2,g即可任取i的因子,最終的g和i/g各乘d即可,所以可以進行如下化簡。(關鍵①)
$$F(n)=\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d^2}}\sum_{g|i}[(g,\frac{i}{g})=1]$$
轉化成φ希望不大,所以直接莫比烏斯反演。
$$F(n)=\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d^2}}\sum_{g|i}\sum_{d‘|g \cap d‘|\frac{i}{g}}\mu(d‘)$$
$$F(n)=\sum_{d=1}^{n}\sum_{d‘=1}^{n}d*\mu(d‘)\sum_{i=1}^{\frac{n}{d^2}}\sum_{g|i}[d‘|g \cap d‘|\frac{i}{g}]$$
這裏和上面一樣,都是要求d‘|g&&d‘|i/g,因此從i中提取2個d‘,即令i=i/d‘^2。
$$F(n)=\sum_{d=1}^{n}\sum_{d‘=1}^{n}d*\mu(d‘)\sum_{i=1}^{\frac{n}{(dd‘)^2}}\sum_{g|i}1$$
會發現後面是約數個數。(關鍵②)
$$F(n)=\sum_{d=1}^{n}\sum_{d‘=1}^{n}d*\mu(d‘)\sum_{i=1}^{\frac{n}{(dd‘)^2}}\sigma_0(i)$$
前面部分發現d*μ(d‘)好像可以卷積到φ,考慮合並dd‘來構造卷積,令d=dd‘。(關鍵③)
$$F(n)=\sum_{d=1}^{\sqrt{n}}\sum_{g|d}g*\mu(\frac{n}{g})\sum_{i=1}^{\frac{n}{(dd‘)^2}}\sigma_0(i)$$
這裏d只枚舉到√n,因為d>√n時後面的Σ=0,沒有貢獻。
然後根據狄利克雷卷積μ*id=φ可以化簡
$$F(n)=\sum_{d=1}^{\sqrt{n}}\varphi(d)\sum_{i=1}^{\frac{n}{(dd‘)^2}}\sigma_0(i)$$
大功告成!
其中,約數個數的前綴和可以進行分塊取值優化,如下
$$\sum_{i=1}^{n}\sigma_0(i)=\sum_{i=1}^{n}\sum_{d|i}1=\sum_{i=1}^{n}\sum_{j=1}^{\frac{n}{i}}1$$
$$\sum_{i=1}^{n}\sigma_0(i)=\sum_{i=1}^{n}\left \lfloor \frac{n}{i} \right \rfloor$$
然後線性篩φ的過程中求解即可。
復雜度分析:
$$\sum_{i=1}^{\sqrt{n}}O(\sqrt{\frac{n}{i^2}})=\sum_{i=1}^{\sqrt{n}}O(\frac{\sqrt{n}}{i})=O(\sqrt{n} ln \sqrt{n})$$
倒數第二步將√n提到最外面,Σ內就是調和數列,和近似為ln n。
#include<cstdio> #include<cmath> #define int long long const int maxn=32000000; int phi[maxn],n,prime[maxn],tot; int solve(int n){ int pos,sum=0; for(int i=1;i<=n;i=pos+1){ pos=n/(n/i); sum+=(pos-i+1)*(n/i); } return sum; } #undef int int main(){ #define int long long scanf("%lld",&n); int ans=1*solve(n),N=(int)sqrt(n)+1;phi[1]=1;//1 for(int i=2;i<=N;i++){ if(!phi[i])phi[prime[++tot]=i]=i-1; for(int j=1;j<=tot&&i*prime[j]<=N;j++){ if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;} phi[i*prime[j]]=phi[i]*(prime[j]-1); } ans+=phi[i]*solve(n/(i*i)); } printf("%lld",ans); return 0; }View Code
【Project Euler】530 GCD of Divisors 莫比烏斯反演