習題:Crash的數字表格(莫比烏斯反演)
阿新 • • 發佈:2021-01-09
題目
思路
貌似可以直接推式子來做
\[\begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)&=\sum_{i=1}^n\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}\\&=\sum_{t=1}^{min(n,m)}\frac{1}{t}\sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)==t]\\&=\sum_{t=1}^{min(n,m)}\frac{1}{t}\sum_{i=1}^{it\le n}\sum_{j=1}^{jt\le m}itjt[gcd(i,j)==1]\\&=\sum_{t=1}^{min(n,m)}\frac{1}{t}\sum_{i=1}^{it\le n}\sum_{j=1}^{jt\le m}itjt\sum_{d|gcd(i,j)}\mu (d)\\&=\sum_{t=1}^{min(n,m)}t\sum_{i=1}^{it\le n}\sum_{j=1}^{jt\le m}ij\sum_{d|gcd(i,j)}\mu (d)\\&=\sum_{t=1}^{min(n,m)}t\sum_{d=1}^{\frac{min(n,m)}{t}}\mu(d)d^2\sum_{i=1}^{itd\le n}\sum_{j=1}^{jtd\le m}ij\\&=\sum_{t=1}^{min(n,m)}t\sum_{d=1}^{\frac{min(n,m)}{t}}\mu(d)d^2(\sum_{i=1}^{itd\le n}i)(\sum_{j=1}^{jtd\le m}j)\\&=\sum_{t=1}^{min(n,m)}t\sum_{d=1}^{\frac{min(n,m)}{t}}\mu(d)d^2(\sum_{i=1}^{\frac{n}{td}}i)(\sum_{j=1}^{\frac{m}{td}}j)\end{aligned} \]注意到這裡雖然是\(\frac{n}{td}\),但是實際上應該是\(\lfloor\frac{\lfloor\frac{n}{t}\rfloor}{d}\rfloor\)
即可以進行兩次數論分塊,具體而言,
\[設g(n,m)=\sum_{d=1}^{min(n,m)}\mu(d)d^2\sum_{i=1}^{\frac{n}{d}}i\sum_{i=1}^{\frac{m}{d}}i\\那麼原式即為\sum_{t=1}^{min(n,m)}t*g(\frac{n}{t},\frac{m}{t}) \]程式碼
#include<iostream> #include<cmath> using namespace std; const int mod=20101009; const int inv2=10050505; int n,m; bool vis[10000005]; int pri[1000005],lenp; long long mu[10000005]; long long f[10000005]; long long ans; void sieve(int n) { mu[1]=1; f[1]=1; for(int i=2;i<=n;i++) { f[i]=(f[i-1]+i)%mod; if(vis[i]==0) { pri[++lenp]=i; mu[i]=-1; } for(int j=1;j<=lenp&&1ll*i*pri[j]<=n;j++) { vis[pri[j]*i]=1; mu[pri[j]*i]=-mu[i]; if(i%pri[j]==0) { mu[i*pri[j]]=0; break; } } mu[i]=mu[i]*i*i%mod; mu[i]=(mu[i]+mu[i-1])%mod; } } long long getsum(long long x) { return x*(x+1)%mod*inv2%mod; } long long solve(int n,int m) { long long temp=0; for(int l=1,r;l<=min(n,m);l=r+1) { r=min(min(n,m),min(n/(n/l),m/(m/l))); temp=(temp+(mu[r]-mu[l-1])%mod*getsum(n/r)%mod*getsum(m/r)%mod)%mod; } temp=(temp%mod+mod)%mod; return temp; } int main() { cin>>n>>m; sieve(min(n,m)); for(int l=1,r;l<=min(n,m);l=r+1) { r=min(min(n,m),min(n/(n/l),m/(m/l))); //cout<<"debug:"<<solve(n/r,m/r)<<'\n'; ans=(ans+((f[r]-f[l-1])%mod+mod)%mod*solve(n/r,m/r)%mod)%mod; } cout<<ans; return 0; } /* n/td n/i */