【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

#include <cctype>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int MAXSIZE=10000020;
int bufpos;
char buf[MAXSIZE];
void init(){
    #ifdef LOCAL
        freopen("3994.txt","r",stdin);
    #endif // LOCAL
    buf[fread(buf,1,MAXSIZE,stdin)]='\0';
    bufpos=0;
}
int readint(){
    int val=0;
    for(;!isdigit(buf[bufpos]);bufpos++);
    for(;isdigit(buf[bufpos]);bufpos++)
        val=val*10+buf[bufpos]-'0';
    return val;
}
const int maxn=1000000;
int prime[maxn+1],cur=0,u[maxn+1],c[maxn+1],d[maxn+1];
int sumu[maxn+1];
long long g[maxn+1];
bool notprime[maxn+1];
void sieve(){
    notprime[1]=true;
    u[1]=d[1]=1; //c(1) is undefined
    for(int i=2;i<=maxn;i++){
        if (!notprime[i]){
            prime[++cur]=i;
            u[i]=-1;
            c[i]=1;
            d[i]=2;
            //continue;
        }
        for(int j=1;j<=cur && i*prime[j]<=maxn;j++){
            int t=i*prime[j];
            notprime[t]=true;
            if (i%prime[j]==0){
                d[t]=d[i]/(c[i]+1)*(c[i]+2);
                c[t]=c[i]+1;
                u[t]=0;
                break;
            }
            u[t]=-u[i];
            d[t]=d[i]*d[prime[j]];
            c[t]=1;
        }
    }
    for(int i=1;i<=maxn;i++)
        sumu[i]=sumu[i-1]+u[i],g[i]=g[i-1]+d[i];
}
int main(){
    init();
    sieve();
    int t=readint();
    while(t--){
        int n=readint(),m=readint(),w=min(n,m),nxt;
        long long ans=0;
        for(int i=1;i<=w;i=nxt+1){
            nxt=min(n/(n/i),m/(m/i));
            ans+=(sumu[nxt]-sumu[i-1])*g[n/i]*g[m/i];
        }
        printf("%lld\n",ans);
    }
}

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

隐藏