莫比乌斯反演泛做1

1.AT5200 [AGC038C] LCMs

Problem

\(\sum_{i=1}^N\sum_{j=i+1}^Nlcm(A_i,A_j)\)

Sol

考虑先把下标变换到\(\sum_{i=1}^N\sum_{j=1}^N\)的形式，如果下标是这样的形式，那么可以看做是一个\(N\times N\)的矩阵，对角线上的元素是\(lcm(A_i,A_i)=A_i\)，上三角和下三角是对称的即\(lcm(A_i,A_j)=lcm(A_j,A_i)\)，因此有

\[\sum_{i=1}^N\sum_{j=i+1}^Nlcm(A_i,A_j)\\ =\frac{\sum_{i=1}^N\sum_{j=1}^Nlcm(A_i,A_j)-\sum_{i=1}^NA_i}{2}\\ \]

接下来只要考虑\(\sum_{i=1}^N\sum_{j=1}^Nlcm(A_i,A_j)\)即可。

\[\sum_{i=1}^N\sum_{j=1}^Nlcm(A_i,A_j)\\ =\sum_{i=1}^N\sum_{j=1}^N\frac{A_i\times A_j}{gcd(A_i,A_j)}\\ =\sum_{i=1}^N\sum_{j=1}^N\sum_{d=max(A)}^N\frac{A_i\times A_j}{d}[gcd(A_i,A_j)==d]\\ =\sum_{d=1}^{max(A)}\frac{1}{d}\sum_{i=1}^N\sum_{j=1}^NA_iA_j[gcd(A_i,A_j)==d]\\ =\sum_{d=1}^{max(A)}\frac{1}{d}\sum_{i=1}^N\sum_{j=1}^NA_iA_j[d|A_i][d|A_j][gcd(\frac{A_i}{d},\frac{A_j}{d})==1]\\ =\sum_{d=1}^{max(A)}\frac{1}{d}\sum_{i=1}^N\sum_{j=1}^NA_iA_j[d|A_i][d|A_j]\sum_{t|gcd(\frac{A_i}{d},\frac{A_j}{d})}\mu(t)\\ =\sum_{d=1}^{max(A)}\sum_{t=1}^{\lfloor \frac{max(A)}{d}\rfloor}\frac{\mu(t)}{d}\sum_{i=1}^N\sum_{j=1}^NA_iA_j[dt|A_i][dt|A_j]\\ =\sum_{d=1}^{max(A)}\sum_{t=1}^{\lfloor \frac{max(A)}{d}\rfloor}\frac{\mu(t)}{d}(\sum_{i=1}^NA_i[dt|A_i])^2\\ =\sum_{d=1}^{max(A)}\sum_{t=1}^{\lfloor \frac{max(A)}{d}\rfloor}\frac{\mu(t)}{d}(\sum_{dt|i}^{max(A)}i\times cnt_i)^2 \]

后面\(\sum_{dt|i}^{max(A)}i\times cnt_i\)可以预处理，复杂度是\(O(nlogn\)的。

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e6+10;
const int mod=998244353;
int primes[N],st[N],cnt;
int mu[N],sum[N];
ll g[N];
int num[N];
void get_primes()
{
    mu[1]=1;
    for(int i=2;i<=N-10;i++)
    {
        if(!st[i]) primes[++cnt]=i,mu[i]=-1;
        for(int j=1;primes[j]*i<=N-10;j++)
        {
            st[primes[j]*i]=1;
            if(i%primes[j]==0) break;
            mu[i*primes[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N-10;i++) sum[i]=sum[i-1]+mu[i];
}
ll power(ll x,int y)
{
    ll res=1;
    while(y)
    {
        if(y&1) res=res*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return res;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    get_primes();
    int n;
    cin>>n;
    vector<int>a(n+1);
    int mx=0;
    for(int i=1;i<=n;i++) cin>>a[i],mx=max(mx,a[i]),num[a[i]]=(num[a[i]]+a[i])%mod;
    for(int i=1;i<=mx;i++)
        for(int j=i;j<=mx;j+=i)
        g[i]=(g[i]+num[j])%mod;
    ll ans=0;
    for(int d=1;d<=mx;d++)
    {
        ll res=0;
        for(int t=1;t<=mx/d;t++)
        res=(res+mu[t]*g[d*t]%mod*g[d*t]%mod+mod)%mod;
        res=res*power(d,mod-2)%mod;
        ans=(ans+res)%mod;
    }
    for(int i=1;i<=n;i++) ans=(ans-a[i]+mod)%mod;
    ans=ans*power(2,mod-2)%mod;
    ans=(ans+mod)%mod;
    cout<<ans<<'\n';
    return 0;
}

2.P6271 [湖北省队互测 2014]一个人的数论

Problem

计算\(\sum_{i=1}^n[gcd(i,n)==1]i^k\)

Sol

\[\sum_{i=1}^n[gcd(i,n)==1]i^k=\sum_{i=1}^ni^k\sum_{d|i,d|n}\mu(d)\\ 令j=\frac{i}{d}，则\\ \sum_{i=1}^ni^k\sum_{d|i,d|n}\mu(d)=\sum_{i=1}^nd^kj^k\sum_{d|i,d|n}\mu(d)=\sum_{d|n}d^k\mu(d)\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}j^k \]

后面是一个自然数幂和的形式，解决办法有第二类斯特林数、伯努利数、拉格朗日插值等，这里考虑使用伯努利数来解决。

• 伯努利数
  • 定义
    • 递归定义：\(B_0=1\)，\(B_n=-\frac{1}{n+1}\sum_{k=0}^{n-1}\binom{n+1}{k}B_k\)
    • 生成函数定义：\(f(x)=\frac{x}{e^x-1}=\sum_{i=0}^{\infty}B_i\frac{x^i}{i!}\)
• 自然数幂和
  \[\sum_{i=1}^ni^k=\frac{1}{1+k}\sum_{i=0}^{k}\binom{k+1}{i}B_in^{k+1-i}=\frac{1}{1+k}\sum_{i=1}^{k+1}\binom{k+1}{i}B_{k+1-i}(n+1)^i\\ 令f_i=\frac{\binom{k+1}{i}B_i}{k+1}，则\sum_{i=1}^ni^k=\sum_{i=0}^{k}f_in^{k+1-i} \]
  于是
  \[\sum_{d|n}d^k\mu(d)\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}j^k \\ =\sum_{d|n}d^k\mu(d)\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor -1}j^k+\sum_{d|n}d^k\mu(d){\lfloor \frac{n}{d} \rfloor}^{k}\\ =\sum_{d|n}d^k\mu(d)\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor -1}j^k+n^k\sum_{d|n}\mu(d)\\ =\sum_{d|n}d^k\mu(d)\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor -1}j^k+n^k[n==1]\\ =\sum_{d|n}d^k\mu(d)\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor -1}j^k\\ =\sum_{d|n}d^k\mu(d)\sum_{i=0}^{k}f_i(\frac{n}{d})^{k+1-i}\\ =\sum_{i=0}^{k}f_in^{k+1-i}\sum_{d|n}d^{i-1}\mu(d)\\ 设h(n)=id(n)_{i-1}\times \mu(n)，g(n)=I(n)=1，其中id(n)_{i-1}=n^{i-1}，则\\ F(n)=h(n)*g(n)=\sum_{d|n}h(d)g(\frac{n}{d})必为积性函数\\ 由于积性函数F(n)=\prod_iF(p_i^{\alpha_i})\\ 对于F(p_i^{\alpha_i})=\sum_{d|p_i^{\alpha_i}}\mu(d)d^{i-1},d=1,p_i,...,p_i^{\alpha_i}，当且仅当d=1和d=p_i时，\mu(d)不为0\\ 于是F(p_i^{\alpha_i})=\mu(1)+\mu(p_i)p_i^{i-1}=1-p_i^{i-1}\\ 综上\quad \sum_{i=0}^{k}f_in^{k+1-i}\sum_{d|n}d^{i-1}\mu(d)=\quad \sum_{i=0}^{k}f_in^{k+1-i}F(n)=\quad \sum_{i=0}^{k}f_in^{k+1-i}\prod_j(1-p_j^{i-1}) \]

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod=1e9+7;
const int N=105;
ll fac[N],infac[N];
ll power(ll x,int y)
{
  ll res=1;
  while(y)
  {
    if(y&1) res=res*x%mod;
    x=x*x%mod;
    y>>=1;
  }
  return res;
}
ll C(int x,int y){
  return fac[x]*infac[y]%mod*infac[x-y]%mod;
}
struct Bernou{
  int d;
  vector<ll>B,f;
  Bernou(int d_):d(d_),B(d+2),f(d+2){}
  void init()
  {
    B[0]=1;
    ll sum=0;
    for(int n=1;n<=d;n++)
    {
        ll sum=0;
        for(int k=0;k<n;k++) sum=(sum+C(n+1,k)*B[k]%mod)%mod;
        B[n]=(-power(n+1,mod-2)*sum%mod+mod)%mod;
    }
    ll inv=power(d+1,mod-2);
    for(int i=0;i<=d;i++) f[i]=C(d+1,i)*B[i]%mod*inv%mod;
  }
};
int main()
{
  ios::sync_with_stdio(false);
  int d,w;
  cin>>d>>w;
  fac[0]=1;
  for(int i=1;i<=d+1;i++) fac[i]=fac[i-1]*i%mod;
  infac[d+1]=power(fac[d+1],mod-2)%mod;
  for(int i=d+1;i>=1;i--) infac[i-1]=infac[i]*i%mod;
  Bernou ber(d);
  ber.init();
  vector<ll>p(w+1),a(w+1);
  ll n=1;
  for(int i=1;i<=w;i++) cin>>p[i]>>a[i],n=n*power(p[i],a[i])%mod;
    //cout<<n<<'\n';
  ll ans=0;
  // for(int i=0;i<=d;i++) cout<<ber.B[i]<<" \n"[i==d];
  // for(int i=0;i<=d;i++) cout<<ber.f[i]<<" \n"[i==d];
  for(int i=0;i<=d;i++)
  {
     ll prod=1;
     ll pw;
     if(i==0) pw=mod-2;
     else pw=i-1;
     for(int j=1;j<=w;j++)prod=prod*(1-power(p[j],pw)+mod)%mod;
    ans=(ans+ber.f[i]*power(n,d+1-i)%mod*prod%mod)%mod;
  }
  cout<<ans<<'\n';
  return 0;
}

3.YY 的 GCD[莫比乌斯反演、\(\mu\)函数求前缀和优化]

Problem

给定\(N\)，\(M\)，求\(1\le x\le N\)，\(1\le y\le M\)且\(gcd(x,y)\)为质数的\((x,y\)有多少对

Sol

\[\sum_{k=1}^{min(n,m)}\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==k][k\in primes]\\ =\sum_{k=1}^{min(n,m)}\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor}[gcd(i,j)==1][k\in primes]\\ =\sum_{k=1}^{min(n,m)}\sum_{i=1}^{\lfloor \frac{n}{k} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{k} \rfloor}\sum_{d=1}^{min(\lfloor \frac{n}{k} \rfloor,\lfloor \frac{m}{k} \rfloor)}\mu(d)[k\in primes][d|i][d|j]\\ =\sum_{k=1}^{min(n,m)}\sum_{d=1}^{min(\lfloor \frac{n}{k} \rfloor,\lfloor \frac{m}{k} \rfloor)}\mu(d){\lfloor \frac{n}{kd} \rfloor} {\lfloor \frac{m}{kd} \rfloor} [k\in primes]\\ =\sum_{T=1}^{min(n,m)}{\lfloor \frac{n}{T} \rfloor} {\lfloor \frac{m}{T} \rfloor}\sum_{k=1}\mu(\frac{T}{k})[k\in primes\ \&\& \ k|T]\ \ \ \ let\ \ \ T=kd\\ \]

后面\(\sum_{k=1}\mu(\frac{T}{k})\)可以用线性筛快速求出关于\(T\)的前缀答案，注意到虽然求的是\(\mu(\frac{T}{k})\)，但对于一个固定的\(T\)，贡献是不变的，所以我们就计算\(T\)的贡献就好了。然后前面用整除分块即可。

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e7+10;
int primes[N],st[N],cnt,mu[N],sum[N],f[N];
void get_primes()
{
    mu[1]=1;
    for(int i=2;i<=10000000;i++)
    {
        if(!st[i]) primes[++cnt]=i,mu[i] = -1;
        for(int j=1;primes[j]*i<=10000000;j++)
        {
            st[primes[j]*i]=1;
            if(i%primes[j]==0) break;
            mu[primes[j] * i] = -mu[i];
        }
    }
    for(int i=1;i<=cnt;i++)
        for(int j=1;j*primes[i]<=10000000;j++)
            f[j*primes[i]]+=mu[j];
    for (int i = 1;i<=10000000;i ++ ) sum[i] = sum[i - 1] + f[i];
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    get_primes();
    int T;
    cin>>T;
    while(T--)
    {
        int n,m;
        cin>>n>>m;
        if(n>m) swap(n,m);
        ll ans=0;
        for(int T=1,nT=0;l<=n;T=nT+1)
        {
            nT=min(n/(n/T),m/(m/T));
            ans+=1LL*(sum[nT]-sum[T-1])*(n/T)*(m/T);
        }
        cout<<ans<<'\n';
    }
}

4.[HAOI2011]Problem b[莫比乌斯反演]

Problem

对于给出的 \(n\) 个询问，每次求有多少个数对 \((x,y)\)，满足 \(a \le x \le b\)，\(c \le y \le d\)，且 \(gcd(x,y) = k\)。

数据范围都在\(5\times 10^4\)

Solve

\[f(n,m,k)=\sum_{i=1}^nsum_{j=1}^ngcd(i,j)=k\\ =\sum_{d=1}^{min(n,m)}\mu(d){\lfloor \frac{n}{kd} \rfloor}{\lfloor \frac{m}{kd} \rfloor} \]

然后由于下标不是从\(1\)开始的，所以考虑容斥：\(res=f(b,d,k)-f(b,c-1,k)-f(a-1,d,k)+f(a-1,c-1,k)\)

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e4+10;
int primes[N],st[N],mu[N],cnt,sum[N];
void get_primes()
{
    mu[1]=1;
    for(int i=2;i<=50000;i++)
    {
        if(!st[i]) primes[++cnt]=i,mu[i]=-1;
        for(int j=1;primes[j]*i<=50000;j++)
        {
            st[primes[j]*i]=1;
            if(i%primes[j]==0) break;
            mu[primes[j]*i]=-mu[i];

        }
    }
    for(int i=1;i<=50000;i++) sum[i]=sum[i-1]+mu[i];
}
ll cal(int n,int m,int k)
{
    ll ans=0;
    if(n>m) swap(n,m);
    for(int l=1,r;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        ans+=1LL*(sum[r]-sum[l-1])*(1LL*n/(1LL*l*k))*(1LL*m/(1LL*l*k));
    }
    return ans;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    get_primes();
    int T;
    cin>>T;
    while(T--)
    {
        int a,b,c,d,k;
        cin>>a>>b>>c>>d>>k;
        ll res=cal(b,d,k)-cal(b,c-1,k)-cal(a-1,d,k)+cal(a-1,c-1,k);
        cout<<res<<'\n';
    }
}

5.P3768-简单的数学题[莫比乌斯反演,杜教筛]

Problem

给定\(n\)，\(p\)，计算\((\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j))mod\ p\)

\(1\le n\le 10^{10}\)，\(5\times 10^8 \le p \le 1.1\times 10^9\)，\(p\)为质数

Solve

\[\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)\\ =\sum_{d}^nd\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)==d]\\ =\sum_{d}^nd^3\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}ij[gcd(i,j)==1]\\ =\sum_{d}^nd^3\sum_{k=1}^{\lfloor \frac{n}{d} \rfloor}\mu(k)\sum_{i=1}^{\lfloor \frac{n}{dk} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{dk} \rfloor}ik\times jk\\ =\sum_{d}^nd^3\sum_{k=1}^{\lfloor \frac{n}{d} \rfloor}\mu(k){k}^2(1+2+...+\lfloor \frac{n}{dk} \rfloor)^2\\ =\sum_{d}^nd^3\sum_{k=1}^{\lfloor \frac{n}{d} \rfloor}\mu(k){k}^2(1+2+...+\lfloor \frac{n}{dk} \rfloor)^2\\ 令T=dk，则\\ 上式=\sum_{T=1}^n(1+2+...+\lfloor \frac{n}{T} \rfloor)^2\sum_{d|T}d^3\mu(\frac{T}{d})(\frac{T}{d})^2\\ =\sum_{T=1}^n(1+2+...+\lfloor \frac{n}{T} \rfloor)^2T^2\sum_{d|T}d\mu(\frac{T}{d}) \]

注意到后面\(T^2\sum_{d|T}d\mu(\frac{T}{d})=T^2\mu*id(T)=T^2\varphi(T)\)

\(\varphi=\epsilon*\varphi=\mu*I*\varphi=\mu*id\)

现在的问题是如何快速求\(\sum_{T=1}^nT^2\varphi(T)\)，用杜教筛快速求

\[S(n)=\sum_{i=1}^ni^2\varphi(i)=\sum_{i=1}^ni^2\varphi(i)= \sum_{i=1}^n(id^2\varphi)*I(i)=\sum_{i=1}^n 1\times S({\lfloor \frac{n}{i} \rfloor})\\ id^2*S(n)=\sum_{i=1}^n(id^2\varphi)*id^2(i)=\sum_{i=1}^nid^2(i)S({\lfloor \frac{n}{i} \rfloor})\\ S(n)=\sum_{i=1}^n(id^2\varphi)*id^2(i)-\sum_{i=2}^nid^2(i)S({\lfloor \frac{n}{i} \rfloor}) \]

现在考虑\((id^2\varphi)*id^2(i)\)

\[(id^2\varphi)*id^2(i)=\sum_{d|i}(id^2\varphi)(d)id^2(\frac{i}{d})\\ =\sum_{d|i}d^2\varphi(d)(\frac{i}{d})^2=\sum_{d|i}\varphi(d)i^2\\ =i^2\sum_{d|i}\varphi(d)=i^3 \]

于是

\[S(n)=\sum_{i=1}^ni^3-\sum_{i=2}^ni^2S({\lfloor \frac{n}{i} \rfloor})\\ =(\frac{n(n+1)}{2})^2-\sum_{i=2}^ni^2S({\lfloor \frac{n}{i} \rfloor})\\ 而\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6} \]

因此最后的答案就是

\[res=\sum_{T=1}^n(\frac{(\lfloor \frac{n}{T} \rfloor)(\lfloor \frac{n}{T} \rfloor+1)}{2})^2T^2\sum_{d|T}d\mu(\frac{T}{d}) \]

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=8e6+10;
map<ll,ll>Sphi;
int primes[N],st[N],cnt;
ll phi[N];
ll inv2,inv6,p;
void get_primes()
{
    phi[1]=1;
    for(int i=2;i<=N-10;i++)
    {
        if(!st[i]) primes[++cnt]=i,phi[i]=i-1;
        for(int j=1;primes[j]*i<=N-10;j++)
        {
            st[primes[j]*i]=1;
            if(i%primes[j]==0)
            {
                phi[i*primes[j]]=phi[i]*primes[j];
                break;
            }
            phi[i*primes[j]]=phi[i]*(primes[j]-1);
        }
    }
    for(int i=1;i<=N-10;i++) phi[i]=(phi[i-1]+phi[i]*i%p*i%p)%p;//预处理部分i^2phi[i]

}
ll s1(ll x){
    x%=p;
   return x*(x+1)%p*inv2%p;
}
ll  s2(ll x){
    x%=p;
 return x*(x+1)%p*(x+x+1)%p*inv6%p;
}
ll S(ll n)
{

    if(n<=N-10) return phi[n];
    if(Sphi[n]) return Sphi[n];
     ll res=s1(n);
     res=res*res%p;
     for(ll l=2,r;l<=n;l=r+1)
     {
        r=n/(n/l);
        res=(res-(s2(r)-s2(l-1)+p)%p*S(n/l)%p+p)%p;
     }
     return Sphi[n]=(res+p)%p;
}
ll power(ll x,ll y,ll P)
{
    ll res=1;
    while(y)
    {
       if(y&1) res=res*x%P;
       x=x*x%P;
       y>>=1;
     }
     return res;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    //freopen("out1.out","w",stdout);
    ll n;
    cin>>p>>n;
    get_primes();
    ll ans=0;
    inv2=power(2,p-2,p);
    inv6=power(6,p-2,p);
    for(ll l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        ll t=s1(n/l);
        t=t*t%p;
        ans=(ans+(S(r)-S(l-1)+p)%p*t%p)%p;
    }
    ans=(ans+p)%p;
    cout<<ans<<'\n';
}

6.[SDOI2015]约数个数和[莫比乌斯反演]

Problem

设\(d(x)\)表示\(x\)的约数个数，给定\(n,m\)，求\(\sum_{i=1}^n\sum_{j=1}^md(ij)\)

Solve

利用\(d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\)

\[\sum_{i=1}^n\sum_{j=1}^md(ij)=\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\\ =\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}\sum_{d=1}^{min(n,m)}\mu(d)[d|gcd(x,y)]\\ =\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[d|gcd(x,y)]\\ =\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^n\sum_{y=1}^m[d|gcd(x,y)]\sum_{i=1}^{\lfloor \frac{n}{x} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{y} \rfloor} \ (变为枚举约数x,y)\\ =\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^n\sum_{y=1}^m[d|gcd(x,y)]{\lfloor \frac{n}{x} \rfloor}{\lfloor \frac{m}{y} \rfloor}\\ =\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{y} \rfloor}{\lfloor \frac{n}{dx} \rfloor}{\lfloor \frac{m}{dy} \rfloor}\\ =\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor \frac{n}{d} \rfloor}{\lfloor \frac{n}{dx} \rfloor}\sum_{y=1}^{\lfloor \frac{m}{y} \rfloor}{\lfloor \frac{m}{dy} \rfloor} \]

后面部分可以预处理出来

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e4+10;
int primes[N],st[N],cnt;
int mu[N],sum[N];
ll g[N];
void get_primes()
{
    mu[1]=1;
    for(int i=2;i<=N-10;i++)
    {
        if(!st[i]) primes[++cnt]=i,mu[i]=-1;
        for(int j=1;primes[j]*i<=N-10;j++)
        {
            st[primes[j]*i]=1;
            if(i%primes[j]==0) break;
            mu[i*primes[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N-10;i++) sum[i]=sum[i-1]+mu[i];

    for(int i=1;i<=N-10;i++)
    {
         ll ans=0;
         for(int l=1,r;l<=i;l=r+1)
         {
            r=i/(i/l);
            ans+=1LL*(r-l+1)*(i/l);
         }
         g[i]=ans;
    }
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    get_primes();
    int T;
    cin>>T;
    while(T--)
    {
         int n,m;
         cin>>n>>m;
         if(n>m) swap(n,m);
         ll ans=0;
         for(int l=1,r;l<=n;l=r+1)
         {
            r=min(n/(n/l),m/(m/l));
            ans+=1LL*(sum[r]-sum[l-1])*g[n/l]*g[m/l];

         }
         cout<<ans<<'\n';
    }
    return 0;
}