- admin 的博客
1517: *【莫比乌斯反演】i*j*gcd(i,j)求和 [scy]
- @ 2025-10-13 11:45:07
问题
给出 ,求 $\sum\limits_{i=1}^n\sum\limits_{j=1}^m i*j*\gcd(i,j)$ 。
题解
| 设 | 111 |
|---|---|
| $\sum\limits_{i=1}^n\sum\limits_{j=1}^m i*j*\gcd(i,j)$ | 222 |
| $=\sum\limits_{i=1}^n\sum\limits_{j=1}^m \sum\limits_{d=1}^n i*j*d \lbrack\gcd(i,j)=d\rbrack$ | 333 |
| $=\sum\limits_{d=1}^n \sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{m}{d}} id*jd*d \lbrack\gcd(i,j)=1\rbrack$ | 444 |
| $=\sum\limits_{d=1}^n d^3 \sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{m}{d}} i*j \lbrack \gcd(i,j)=1 \rbrack$ | 555 |
| $ =\sum\limits_{d=1}^n d^3 \sum\limits_{i=1}^ {\frac_n_d} \sum\limits_{j=1}^{\frac{m}{d}}i*j \sum \limits_{a | \gcd(i,j)} \ mu(a) $ |
| $=\sum\limits_{d=1}^n d^3 \sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{m}{d}}i*j\sum\limits_{a=1}^{\frac{n}{d}} \lbrack a | i \rbrack \lbrack a |
| $=\sum\limits_{d=1}^n d^3 \sum\limits_{a=1}^{\frac{n}{d}}\mu(a)\sum\limits_{i=1}^{\frac{n}{d}}i*\lbrack a | i \rbrack\sum\limits_{j=1}^{\frac{m}{d}}j*\lbrack a |
| $=\sum\limits_{d=1}^n d^3 \sum\limits_{a=1}^{\frac{n}{d}}a^2\mu(a)\sum\limits_{i=1}^{\frac{n}{ad}}i\sum\limits_{j=1}^{\frac{m}{ad}}j$ | 设: |
| $=\sum\limits_{d=1}^n d^3 \sum\limits_{a=1}^{\frac{n}{d}}a^2\mu(a)S(\lfloor \frac{n}{ad} \rfloor) S(\lfloor \frac{m}{ad} \rfloor)$ | 设: |
| $=\sum\limits_{d=1}^n d \sum\limits_{\frac{T}{d}=1}^{\frac{n}{d}} T^2 \mu(\frac{T}{d})S(\lfloor \frac{n}{T} \rfloor)S(\lfloor \frac{m}{T} \rfloor)$ | 777 |
| $=\sum\limits_{d=1}^n d \sum\limits_{T=1}^n T^2 \mu(\frac{T}{d}) S(\lfloor \frac{n}{T} \rfloor)S(\lfloor \frac{m}{T} \rfloor)$ | 888 |
| $= \sum\limits_{T=1}^n T^2 \sum\limits_{d=1}^n d \mu(\frac{T}{d}) S(\lfloor \frac{n}{T} \rfloor)S(\lfloor \frac{m}{T} \rfloor)$ | 设: $=\sum\limits_{d |
| $=\sum\limits_{T=1}^n S(\lfloor \frac{n}{T} \rfloor)S(\lfloor \frac{m}{T} \rfloor) T^2 F(T)$ | 到此,问题解决。但代码运行时间2.23s(时限2.5s,刚好ac),肯定还有优化,也是一个有难度且非常重要的优化。 |
| 观察 $F(T)=\sum\limits_{d | T}d* \mu(\frac{T}{d})$ |
具体代码:
2.3s代码如下
#include<bits/stdc++.h>//1.86s
#define LL long long
using namespace std;
const int N=1e7;
const LL mod=20101009;
int pr, p[N+10];LL mu[N+10],F[N+10],S[N+10];bool v[N+10];
void init()
{
memset(v, 0, sizeof(v));
pr=0;mu[0]=0;mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!v[i]) p[++pr]=i,mu[i]=-1;
for(int j=1;j<=pr&&p[j]*i<=N;j++)
{
v[i*p[j]]=true;
if(i%p[j]==0) {mu[i*p[j]]=0; break;}
mu[i*p[j]]=-mu[i];
}
}
memset(F,0,sizeof(F));
for(int i=1;i<=N;i++)for(int j=i;j<=N;j+=i)
F[j]=(F[j]+mu[j/i]*i%mod+mod)%mod;
for(int i=2;i<=N;i++)F[i]=(F[i]*i%mod*i%mod+F[i-1]+mod)%mod;
S[1]=1;for(int i=2;i<=N;i++)S[i]=(S[i-1]+i)%mod;
}
LL calc(LL n,LL m)
{
LL ans=0;
for(LL l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans=(ans+(F[r]-F[l-1]+mod)%mod*S[n/l]%mod*S[m/l]%mod)%mod;
}
return ans;
}
int main()
{
init();
LL n,m;scanf("%lld%lld",&n,&m);if(n>m)swap(n,m);
printf("%lld\n",calc(n,m));
return 0;
}
0.4s代码如下:
#include<bits/stdc++.h>//400ms
#define LL long long
using namespace std;
const int N=1e7;
const LL mod=20101009;
int pr, p[N+10];LL mu[N+10],F[N+10],S[N+10];bool v[N+10];
void init()
{
memset(v, 0, sizeof(v));memset(F,0,sizeof(F));
pr=0;mu[0]=0;mu[1]=1,F[1]=1;
for(int i=2;i<=N;i++)
{
if(!v[i]) p[++pr]=i,mu[i]=-1,F[i]=(1*mu[i]+i*mu[1])%mod;
for(int j=1;j<=pr&&p[j]*i<=N;j++)
{
v[i*p[j]]=true;
if(i%p[j]==0) {mu[i*p[j]]=0;F[i*p[j]]=F[i]*p[j]%mod; break;}
mu[i*p[j]]=-mu[i];
F[i*p[j]]=F[i]*F[p[j]]%mod;
}
}
for(int i=2;i<=N;i++)F[i]=(F[i]*i%mod*i%mod+F[i-1]+mod)%mod;
S[1]=1;for(int i=2;i<=N;i++)S[i]=(S[i-1]+i)%mod;
}
LL calc(LL n,LL m)
{
LL ans=0;
for(LL l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans=(ans+(F[r]-F[l-1]+mod)%mod*S[n/l]%mod*S[m/l]%mod)%mod;
}
return ans;
}
int main()
{
init();
LL n,m;scanf("%lld%lld",&n,&m);if(n>m)swap(n,m);
printf("%lld\n",calc(n,m));
return 0;
}