P1829 [国家集训队]Crash的数字表格 / JZPTAB
斯德哥尔摩
2018-07-18 00:23:42
[P1829 [国家集训队]Crash的数字表格 / JZPTAB](https://www.luogu.org/problemnew/show/P1829)
题目要求：$Ans=\sum_{i=1}^n\sum_{j=1}^mlcm(i,j) \mod p$
跟某题好像啊：[P3768 简单的数学题](https://www.luogu.org/blog/49998/p3768-jian-dan-di-shuo-xue-ti)
把那个$\mod p$先丢到一边去。
$Ans=\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)$
不妨设：$n<=m$
我们知道：$lcm(i,j)=\frac{ij}{gcd(i,j)}$
带入得：$Ans=\sum_{i=1}^n\sum_{j=1}^m\frac{ij}{gcd(i,j)}$
$\Rightarrow Ans=\sum_{i=1}^n\sum_{j=1}^m\sum_{d|n}\frac{ij}{d}[gcd(i,j)==d]$
$\Rightarrow Ans=\sum_{d=1}^n\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}\frac{idjd}{d}[gcd(i,j)==1]$
能约的都约掉：
$Ans=\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}ij[gcd(i,j)==1]$
把$[gcd(i,j)==1]$反演一下：
$Ans=\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}\sum_{d'|d}id'jd'\mu(d')$
$\Rightarrow Ans=\sum_{d=1}^nd\sum_{d'=1}^{\lfloor \frac{n}{d}\rfloor}\mu(d')(d')^2\sum_{i=1}^{\lfloor \frac{n}{dd'}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{dd'}\rfloor}ij$
在洛谷上，我们的推式子到此结束，并不需要再往后推了。
~~因为开了O2优化~~
然后，线性筛$O(n)$筛出$\sum_{d'=1}^{\lfloor \frac{n}{d}\rfloor}\mu(d')(d')^2$。
后面的那个式子：
$\sum_{i=1}^{\lfloor \frac{n}{dd'}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{dd'}\rfloor}ij=\sum_{i=1}^{\lfloor \frac{n}{dd'}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{dd'}\rfloor}j$
直接等差数列求和即可。
附代码：
```cpp
#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 10000010
#define MOD 20101009LL
using namespace std;
int n,m;
int k=0,prime[MAXN],mu[MAXN];
long long inv2,sum[MAXN];
bool np[MAXN];
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
inline long long SUM(long long x){x%=MOD;return x*(x+1)%MOD*inv2%MOD;}
void make(){
int m=MAXN-10;
mu[1]=1;
for(int i=2;i<=m;i++){
if(!np[i]){
prime[++k]=i;
mu[i]=-1;
}
for(int j=1;j<=k&&prime[j]*i<=m;j++){
np[prime[j]*i]=true;
if(i%prime[j]==0)break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=m;i++)sum[i]=(sum[i-1]+mu[i]*i%MOD*i%MOD)%MOD;
}
long long solve(long long n,long long m){
if(n>m)swap(n,m);
long long s=0;
for(long long i=1,last=1;i<=n;i=last+1){
last=min(n/(n/i),m/(m/i));
s=(s+(sum[last]-sum[i-1])%MOD*SUM(n/i)%MOD*SUM(m/i)%MOD)%MOD;
}
return s;
}
void work(){
long long ans=0;
for(int i=1;i<=n;i++){
ans=(ans+solve(n/i,m/i)*i%MOD)%MOD;
}
printf("%lld\n",(ans+MOD)%MOD);
}
void init(){
n=read();m=read();
inv2=(long long)(MOD+1)/2;
if(n>m)swap(n,m);
}
int main(){
make();
init();
work();
return 0;
}
```
但是，我们如果要达到正解的水平，应该继续推：
设$D=dd'$
$\Rightarrow Ans=\sum_{D=1}^n\sum_{d|D}\mu(\frac{D}{d})(\frac{D}{d})^2d\sum_{i=1}^{\lfloor \frac{n}{D}\rfloor}\sum_{i=1}^{\lfloor \frac{m}{D}\rfloor}ij$
设$sum(n)=1+2+3+...+n$
$\Rightarrow Ans=\sum_{D=1}^nsum(\lfloor \frac{n}{D}\rfloor)sum(\lfloor \frac{m}{D}\rfloor)\sum_{d|D}\mu(\frac{D}{d})\frac{D^2}{d}$
又因为：$d'=\frac{D}{d}$
带入得：$Ans=\sum_{D=1}^nsum(\lfloor \frac{n}{D}\rfloor)sum(\lfloor \frac{m}{D}\rfloor)D\sum_{d'|D}\mu(d')d'$
把$d'$换个名字$x$得：
$Ans=\sum_{D=1}^nsum(\lfloor \frac{n}{D}\rfloor)sum(\lfloor \frac{m}{D}\rfloor)D\sum_{x|D}x\mu(x)$
前面这玩意直接暴力求，关键是后面那一堆$\sum_{x|D}x\mu(x)$
设$f(n)=\sum_{x|n}x\mu(x)$，这是一个积性函数，可以线性筛。
证明如下：
考虑莫比乌斯反演的线性筛，我们为一个数$n$加入一个新的质因子$p$：
如果$n$有这个质因子的话，那么对于$n$原来的约数$x$，这个数如果有质因子$p$，加上这个$p$之后，并不会产生新的质因子。
否则，$μ(px)=0$，所以此时$f(np)=f(n)$。
如果$n$没有这个因子，那么其每一个约数乘上这个因子后，由于$μ(px)=-μ(x)$，所以$f(np)=f(n)-pf(n)=(1-p)f(n)=f(p)f(n)$
附代码：
```cpp#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 10000010
#define MOD 20101009LL
using namespace std;
int n,m;
int k=0,prime[MAXN];
long long s[MAXN],f[MAXN];
bool np[MAXN];
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
void make(){
f[1]=1;
for(int i=2;i<=m;i++){
if(!np[i]){
prime[++k]=i;
f[i]=(MOD+1-i)%MOD;
}
for(int j=1;j<=k&&prime[j]*i<=m;j++){
np[prime[j]*i]=true;
if(i%prime[j]==0){
f[prime[j]*i]=f[i];
break;
}
else f[prime[j]*i]=f[i]*f[prime[j]]%MOD;
}
}
for(int i=1;i<=m;i++){
f[i]=(f[i-1]+f[i]*i%MOD)%MOD;
s[i]=(s[i-1]+i%MOD)%MOD;
}
}
int main(){
long long ans=0;
n=read();m=read();
if(n>m)swap(n,m);
make();
for(int i=1,last=1;i<=n;i=last+1){
last=min(n/(n/i),m/(m/i));
ans=(ans+(f[last]-f[i-1]+MOD)%MOD*s[n/i]%MOD*s[m/i]%MOD)%MOD;
}
printf("%lld\n",ans);
return 0;
}
```