【BZOJ3994】【Vijos1949】【SDOI2015】约数个数和

Author Avatar
q234rty 8月 03, 2016

莫比乌斯反演即可。

由于 NM 的约数一定可以写成 \frac{pM} {q}(p \mid N ,q \mid M) 的形式,我们有：

d(NM)=\sum_{p\mid N}\sum_{q\mid M}[\gcd(p,q)=1]

代入原式：

\begin{aligned}\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)&=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{p \mid i}\sum_{q \mid j}[\gcd(p,q)=1] \\ &=\sum_{i=1}^{N}\sum_{j=1}^{M}\left \lfloor \frac {N} {i} \right \rfloor \left \lfloor \frac {M} {j} \right \rfloor[\gcd(i,j)=1] \\ &=\sum_{i=1}^{N}\sum_{j=1}^{M}\left \lfloor \frac {N} {i} \right \rfloor \left \lfloor \frac {M} {j} \right \rfloor\sum_{d \mid \gcd(i,j)}\mu(d) \\ &=\sum_{d=1}^{N} \mu(d) \sum_{i=1}^{\left \lfloor \frac {N} {d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac {M} {d} \right \rfloor }\left \lfloor\frac{N}{id}\right\rfloor\left \lfloor\frac{M}{jd}\right\rfloor \\ &=\sum_{d=1}^{N} \mu(d) \sum_{i=1}^{\left \lfloor \frac {N} {d} \right \rfloor} \left \lfloor\frac{N}{id}\right\rfloor\sum_{j=1}^{\left \lfloor \frac {M} {d} \right \rfloor }\left \lfloor\frac{M}{jd}\right\rfloor\end{aligned}

设 g(n)=\sum_{i=1}^{n}\left \lfloor \frac {n} {i}\right \rfloor ，于是

\begin{aligned}\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)&=\sum_{d=1}^{N} \mu(d) \sum_{i=1}^{\left \lfloor \frac {N} {d} \right \rfloor} \left \lfloor\frac{\left \lfloor\frac{N}{d}\right\rfloor}{i}\right\rfloor \sum_{j=1}^{\left \lfloor \frac {M} {d} \right \rfloor} \left \lfloor\frac{\left \lfloor\frac{M}{d}\right\rfloor}{j}\right\rfloor\\&=\sum_{d=1}^{N}\mu(d)g(\left \lfloor\frac{N}{d}\right \rfloor) g(\left \lfloor\frac{M}{d}\right\rfloor)\end{aligned}

注意到 g 是 d 的前缀和，而 d 是积性函数，可以线性筛出 d 之后预处理出 g 。

之后用分块加速计算就可以了。

Code

  1. #include <cctype>
  2. #include <cstdio>
  3. #include <cstring>
  4. #include <algorithm>
  5. using namespace std;
  6. const int MAXSIZE=10000020;
  7. int bufpos;
  8. char buf[MAXSIZE];
  9. void init(){
  10.     #ifdef LOCAL
  11.         freopen("3994.txt","r",stdin);
  12.     #endif // LOCAL
  13.     buf[fread(buf,1,MAXSIZE,stdin)]='\0';
  14.     bufpos=0;
  15. }
  16. int readint(){
  17.     int val=0;
  18.     for(;!isdigit(buf[bufpos]);bufpos++);
  19.     for(;isdigit(buf[bufpos]);bufpos++)
  20.         val=val*10+buf[bufpos]-'0';
  21.     return val;
  22. }
  23. const int maxn=1000000;
  24. int prime[maxn+1],cur=0,u[maxn+1],c[maxn+1],d[maxn+1];
  25. int sumu[maxn+1];
  26. long long g[maxn+1];
  27. bool notprime[maxn+1];
  28. void sieve(){
  29.     notprime[1]=true;
  30.     u[1]=d[1]=1; //c(1) is undefined
  31.     for(int i=2;i<=maxn;i++){
  32.         if (!notprime[i]){
  33.             prime[++cur]=i;
  34.             u[i]=-1;
  35.             c[i]=1;
  36.             d[i]=2;
  37.             //continue;
  38.         }
  39.         for(int j=1;j<=cur && i*prime[j]<=maxn;j++){
  40.             int t=i*prime[j];
  41.             notprime[t]=true;
  42.             if (i%prime[j]==0){
  43.                 d[t]=d[i]/(c[i]+1)*(c[i]+2);
  44.                 c[t]=c[i]+1;
  45.                 u[t]=0;
  46.                 break;
  47.             }
  48.             u[t]=-u[i];
  49.             d[t]=d[i]*d[prime[j]];
  50.             c[t]=1;
  51.         }
  52.     }
  53.     for(int i=1;i<=maxn;i++)
  54.         sumu[i]=sumu[i-1]+u[i],g[i]=g[i-1]+d[i];
  55. }
  56. int main(){
  57.     init();
  58.     sieve();
  59.     int t=readint();
  60.     while(t--){
  61.         int n=readint(),m=readint(),w=min(n,m),nxt;
  62.         long long ans=0;
  63.         for(int i=1;i<=w;i=nxt+1){
  64.             nxt=min(n/(n/i),m/(m/i));
  65.             ans+=(sumu[nxt]-sumu[i-1])*g[n/i]*g[m/i];
  66.         }
  67.         printf("%lld\n",ans);
  68.     }
  69. }

知识共享许可协议
本作品采用知识共享署名-相同方式共享 4.0 国际许可协议进行许可。
本文链接：https://www.q234rty.top/2016/08/03/bzoj3994/

隐藏