[luoguP1829][国家集训队]Crash的数字表格 / JZPTAB
求$$sum_{i=1}nsum_{j=1}mLCM(i, j); mod; 20101009$$,
(n, mle {10}^7)
以下(n le m)
[ans=sum_{i=1}^nsum_{j=1}^mLCM(i, j)
]
[=sum_{i=1}^nsum_{j=1}^mfrac{ij}{GCD(i, j)}
]
[=sum_{d=1}^nsum_{i=1}^nsum_{j=1}^mfrac{ij[GCD(i, j) == d]}{d}
]
[=sum_{d=1}^ndsum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{m}{d}
floor}ij[GCD(i, j) == 1]
]
令$$f(d)=sum_{i=1}xsum_{j=1}yij[GCD(i,j) == d]$$
[g(d)=sum_{d|n}f(n)
]
则$$g(d)=sum_{i=1}xsum_{j=1}yij[d|GCD(i,j)]$$
[=d^2sum_{i=1}^{lfloorfrac{x}{d}
floor}sum_{j=1}^{lfloorfrac{y}{d}
floor}ij[1|GCD(i,j)]
]
令(sum(x)=sum_{i=1}^xi)
则$$g(d)=d2sum_{i=1}{lfloorfrac{x}{d}
floor}i*sum(lfloorfrac{y}{d}
floor)$$
[=d^2*sum(lfloorfrac{x}{d}
floor)*sum(lfloorfrac{y}{d}
floor)
]
则$$f(n)=sum_{n|d}mu(frac{d}{n})g(d)$$
[f(1)=sum_{d=1}^xmu(d)g(d)
]
[=sum_{d=1}^xmu(d)d^2*sum(lfloorfrac{x}{d}
floor)*sum(lfloorfrac{y}{d}
floor)
]
[ans=sum_{d=1}^ndsum_{x=1}^{lfloorfrac{n}{d}
floor}mu(x)x^2*sum(lfloorfrac{n}{xd}
floor)*sum(lfloorfrac{m}{xd}
floor)
]
好像到这样就可以了？
void init(){
miu[1]=1;
for(int i=2; i < Maxn; i++){
if(!pr[i]) pr[++ptot]=i, miu[i]=-1;
for(int j=1, x; j <= ptot && (x=i*pr[j]) < Maxn; j++){
pr[x]=1; if(i%pr[j] == 0) break; miu[x]=-miu[i];
}
}
for(ll i=1; i < Maxn; i++) miu[i]=(miu[i-1]+miu[i]*i*i%p)%p, sum[i]=(sum[i-1]+i)%p;
}
ll cal(ll n, ll m){ if(n > m) swap(n, m);
ll ans=0;
for(ll l=1, r=0; r < n; l=r+1){
r=min(n/(n/l), m/(m/l));
ans=(ans+(miu[r]-miu[l-1])*sum[n/l]%p*sum[m/l]%p)%p;
}
return ans;
}
void solve(){
init(); n=read(), m=read(); if(n > m) swap(n, m); ll ans=0;
for(ll l=1, r=0; r < n; l=r+1){
r=min(n/(n/l), m/(m/l));
ans=(ans+1ll*(r+l)*(r-l+1)/2%p*cal(n/l, m/l)%p)%p;
}
printf("%lld
", (ans+p)%p);
}