Luogu P3327 [SDOI2015]约数个数和/BZOJ 4176 Lucas的数论

题目大意

共有 $T$ 组数据，给定 $n,m$，设 $d(x)={\sum }_{i|x}1$，求 $$\sum _{i=1}^n\sum _{j=1}^{m}d(ij)$$

数据范围　$T=1(不输入),1⩽n=m(不输入)⩽10^9$　或　 $1⩽n,m,T⩽50000$

题解

由于两道题目有一定不同，这里考虑最大范围，但不管如何，先来证证证……

---

结论　$\begin{array}{r}d(ij)=\sum _{x|i}^{{}^{{}^{{}^{{}^{}}}}}\sum _{y|j}[(x,y)=1]\end{array}$

证明　对于 $ij$，其唯一分解为 $\begin{array}{r}ij=\prod _{d=1}^{n^{{}^{{}^{{}^{}}}}}{p}_d^{k_d}\end{array}$，则 $\begin{array}{r}d(ij)=\prod _{d=1}^{n^{{}^{{}^{{}^{}}}}}(k_d+1)\end{array}$。
考虑一个质数 $p$ 对 $d(ij)$ 的贡献，设 $i=p^k\ast t_1,j=p^q\ast t_2$，则对 $d(ij)$ 的贡献为 $k+q+1$。
再来看一个引理。

---

引理　设 $i=p^k\ast t_1,j=p^q\ast t_2$，则有：$\begin{array}{r}\sum _{x=0}^k\sum _{y=0}^q[(p^x,p^y)=1]=q+k+1\end{array}$
证明　可以发现，当且仅当 $xy=0$ 时，$(p^x,p^y)=1$，故成立。

---

因此对质数 $p$ 对 $d(i,j)$ 的贡献的贡献为 $\begin{array}{r}\sum _{x=0}^k\sum _{y=0}^q[(p^x,p^y)=1]\end{array}$ 。
接下来设 $\begin{array}{r}i=\prod _{i=1}^n{p}_i^{k_i},j=\prod _{i=1}^{n^{{}^{{}^{{}^{}}}}}{p}_i^{q_i}\end{array}$，考虑下面逐一每一个质数 $p_i$，根据乘法原理得：
$$d(ij)=\sum _{x_1=0}^{k_1}\sum _{y_1=0}^{q_1}[(p_1^{x_1},p_1^{y_1})=1]\sum _{x_2=0}^{k_2}\sum _{y_2=0}^{q_2}[(p_2^{x_2},p_2^{y_2})=1]...\sum _{x_n=0}^{k_n}\sum _{y_n=0}^{q_n}[(p_n^{x_n},p_n^{y_n})=1]$$

即$$d(ij)=\sum _{d=1}^n\sum _{x_d=0}^{k_d}\sum _{y_d=0}^{q_d}[(p_d^{x_d},p_d^{y_d})=1]$$

可以发现，对于下面这个式子
$$\sum _{x|i}\sum _{y|j}[(x,y)=1]$$

如果 $x$ 和 $y$ 含有的不相同的质因子，它们不会对答案造成任何影响，因而相当于仅仅枚举了两者相同的质因子的个数，因此，我们有
$$\sum _{x|i}\sum _{y|j}[(x,y)=1]=\sum _{d=1}^n\sum _{x_d=0}^{k_d}\sum _{y_d=0}^{q_d}[(p_d^{x_d},p_d^{y_d})=1]$$

因而$$d(ij)=\sum _{x|i}\sum _{y|j}[(x,y)=1]$$得证。

---

代入题目，得：
$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^{m}d(ij)& =\sum _{i=1}^n\sum _{j=1}^m\sum _{x|i}\sum _{y|j}[(x,y)=1]\\ & =\sum _{i=1}^n\sum _{j=1}^m\sum _{x|i}\sum _{y|j}\sum _{d|(x,y)}\mu (d)\end{array}$$

各回各家，各找各妈，好好理解，实在不行用C++证明法。
$$\begin{array}{rl}原式& =\sum _{d=1}^{\min (n,m)}\mu (d)\left(\sum _{d|x}^{x⩽n}\sum _{x|i}^{i⩽n}1\right)\left(\sum _{d|y}^{y⩽m}\sum _{y|j}^{i⩽m}1\right)\\ & =\sum _{d=1}^{\min (n,m)}\mu (d)\left(\sum _{d|x}^{x⩽n}\lfloor \frac{n}{x}\rfloor \right)\left(\sum _{d|y}^{y⩽m}\lfloor \frac{m}{y}\rfloor \right)\\ & =\sum _{d=1}^{\min (n,m)}\mu (d)\left(\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{n}{dx}\rfloor \right)\left(\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }\lfloor \frac{m}{dy}\rfloor \right)\\ & =\sum _{d=1}^{\min (n,m)}\mu (d)\left(\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{\lfloor \frac{n}{d}\rfloor }{x}\rfloor \right)\left(\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }\lfloor \frac{\lfloor \frac{m}{d}\rfloor }{y}\rfloor \right)\end{array}$$

Luogu

$\Theta (n\sqrt{n})$预处理出 $f(x)={\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor$，线筛出 $\mu$ 即可，胜在代码短。

#include<bits/stdc++.h>
using namespace std;
#define maxn 50010
long long Sqrt[maxn];
long long s_sqrt(long long n){
    long long ans=0;
    for(long long l=1,r;l<=n;l=r+1){
        r=min(n/(n/l),n);
        ans+=(r-l+1)*(n/l);
    } return ans;
}
long long p[maxn],vis[maxn],mu[maxn],smu[maxn],calced=-1;
void euler(long long n=maxn-10){
    calced=n;
    long long tot=0;
    memset(vis,0,sizeof vis);
    mu[1]=1;
    for(long long i=2;i<=n;i++){
        if(vis[i]==0) p[++tot]=i,mu[i]=-1;
        for(long long j=1;j<=tot&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[p[j]*i]=0;
                break;
            } else mu[p[j]*i]=-mu[i];
        }
    } smu[0]=0;
    for(long long i=1;i<=n;i++)
        smu[i]=smu[i-1]+mu[i];
}
long long T,n,m;
signed main(void){
    euler(50000);
    for(int i=1;i<=maxn;i++)
        Sqrt[i]=s_sqrt(i);
    scanf("%lld",&T);
    while(T--){
        scanf("%lld%lld",&n,&m);
        long long ans=0;
        for(long long l=1,r;l<=min(n,m);l=r+1){
            r=min(min(n/(n/l),m/(m/l)),min(n,m));
            ans+=(smu[r]-smu[l-1])*Sqrt[n/l]*Sqrt[m/l];
        } printf("%lld\n",ans);
    }
    return 0;
}

还能更好吗？这个时间几乎时卡着过去的。当然可以。有
$$\begin{array}{rl}\sum _{i=1}^n\lfloor \frac{n}{i}\rfloor & =\sum _{i=1}^n\sum _{i|k}^{k⩽n}1\\ & =\sum _{k=1}^n\sum _{i|k}1\\ & =\sum _{k=1}^{n}d(k)\end{array}$$

其中 $d$ 是约数和函数。$d$ 可以在 $\Theta (n)$ 的时间内筛出来，时间复杂度为 $\Theta (n+T\sqrt{n})$。

#include<bits/stdc++.h>
using namespace std;
#define maxn 50010
long long low[maxn],d[maxn],sd[maxn];
long long p[maxn],vis[maxn],mu[maxn],smu[maxn],calced=-1;
void euler(long long n=maxn-10){
    calced=n;
    long long tot=0;
    memset(vis,0,sizeof vis);
    mu[1]=1;smu[1]=1;d[1]=1;low[1]=1;sd[1]=1;
    for(long long i=2;i<=n;i++){
        if(vis[i]==0) low[i]=p[++tot]=i,mu[i]=-1,d[i]=2;
        for(long long j=1;j<=tot&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
            	low[i*p[j]]=low[i]*p[j];
                if(low[i]==i) d[i*p[j]]=d[i]+1;
                else d[i*p[j]]=d[i/low[i]]*d[low[i]*p[j]];
                mu[p[j]*i]=0;
                break;
            } mu[p[j]*i]=-mu[i];
			low[i*p[j]]=p[j];
            d[i*p[j]]=d[i]*d[p[j]];
        } smu[i]=smu[i-1]+mu[i];
        sd[i]=sd[i-1]+d[i];
    }
}
long long T,n,m;
signed main(void){
    euler();
    scanf("%lld",&T);
    while(T--){
        scanf("%lld%lld",&n,&m);
        long long ans=0;
        for(long long l=1,r;l<=min(n,m);l=r+1){
            r=min(min(n/(n/l),m/(m/l)),min(n,m));
            ans+=(smu[r]-smu[l-1])*sd[n/l]*sd[m/l];
        } printf("%lld\n",ans);
    }
    return 0;
}

BZOJ

% 　　BZOJ 要求设计非线性算法，因而 $\mu$ 的前缀和需要用杜教筛求出，由于 $n=m$，有：
$$\sum _{d=1}^{\min (n,m)}\mu (d)(\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{\lfloor \frac{n}{d}\rfloor }{x}\rfloor )(\sum _{y=1}^{\lfloor \frac{m}{d}\rfloor }\lfloor \frac{\lfloor \frac{m}{d}\rfloor }{y}\rfloor )=\sum _{d=1}^{\min (n,m)}\mu (d)(\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\lfloor \frac{\lfloor \frac{n}{d}\rfloor }{x}\rfloor {)}^2$$　　这样求第二个sigma的时间复杂度为：$$T(n)=\Theta \left(\sqrt{\lfloor \frac{n}{d}\rfloor }\right)$$　　因而计算这个sigma的总时间复杂度为
$$T(n)=\Theta \left(\sqrt{\frac{n}{1}}\right)+\Theta \left(\sqrt{\frac{n}{2}}\right)+...+\Theta \left(\sqrt{\frac{n}{\sqrt{n}}}\right)=\Theta \left(n^{\frac{3}{4}}\right)$$　　因而总时间复杂度为 $$T(n)=\Theta \left(n^{\frac{2}{3}}+n^{\frac{3}{4}}\right)$$

#include<bits/stdc++.h>
using namespace std;
#define maxn 2000000
#define mod 1000000007
long long s_sqrt(long long n){
    long long ans=0;
    for(long long l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans+(r-l+1)*(n/l))%mod;
    } return ans*ans%mod;
}
long long p[maxn],vis[maxn],mu[maxn],smu[maxn],calced=-1;
void euler(long long n=maxn-10){
    calced=n;
    long long tot=0;
    memset(vis,0,sizeof vis);
    mu[1]=1;smu[1]=1;
    for(long long i=2;i<=n;i++){
        if(vis[i]==0) p[++tot]=i,mu[i]=-1;
        for(long long j=1;j<=tot&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[p[j]*i]=0;
                break;
            } else mu[p[j]*i]=mod-mu[i];
        } smu[i]=(smu[i-1]+mu[i]+mod)%mod;
    }
}
map<long long,long long> Smu;
long long Sum_mu(long long n){
    if(n<=calced) return smu[n];
    map<long long,long long>::iterator it=Smu.find(n);
    if(it!=Smu.end()) return it->second;
    long long ans=1;
    for(long long l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans-Sum_mu(n/l)*(r-l+1)%mod+mod)%mod;
    } return Smu[n]=ans%mod;
}
long long n;
signed main(void){
    euler();
    scanf("%lld",&n);
    long long ans=0;
    for(long long l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans+(Sum_mu(r)-Sum_mu(l-1)+mod)%mod*s_sqrt(n/l)%mod)%mod;
    } printf("%lld\n",ans%mod);
    return 0;
}