Problem

bzoj3994 洛谷3327

題意：設 \(d(x)\) 為 \(x\) 的約數個數，給定 \(N,M\)，求\(\sum_{i=1}^N\sum_{j=1}^Md(ij)\)

\(1\leq N,M,T\leq 5\times 10^4\)

Solution

第一次推出莫反式 ?(^?^*)

以下部分中小括號代表\(\gcd\)，中括號代表取布爾值

一開始想枚舉約數然後看有多少倍數出現過的，發現不好弄，轉而想到 \(xy\) 的因數一定是 \(x\) 的因數和 \(y\) 的因數的積，去除冗余後可以得到 \(d(xy)=\sum_{a|x}\sum_{b|y}[(\frac xa,\frac yb)=1]=\sum_{a|x}\sum_{b|y}[(a,b)=1]\)

然後答案為：

\[\sum_{i=1}^N\sum_{j=1}^M\sum_{a|i}\sum_{b|j}[(a,b)=1]\]

考慮優先枚舉\(a,b\)

\[\sum_{a=1}^N\sum_{b=1}^M[(a,b)=1]\lfloor \frac Na\rfloor \lfloor \frac Mb\rfloor\]

發現裏頭那個\(\gcd\)可以套路地搞事了：

構造\(f(x)=\sum_{a=1}^N\sum_{b=1}^M[(a,b)=x]\lfloor \frac Na\rfloor \lfloor \frac Mb\rfloor\)

構造\(F(x)=\sum_{a=1}^N\sum_{b=1}^M[x|(a,b)]\lfloor \frac Na\rfloor \lfloor \frac Mb\rfloor\)

可以根據定義式得到\(F(x)=\sum_{x|n}f(n)\)，經過莫比烏斯反演，得到\(f(x)=\sum_{x|n}\mu(\frac nx)F(n)\)

答案為\(f(1)=\sum_n\mu(n)F(n)\)

由於\(F(n)\)中可以將\(\gcd\)的限制條件放在前面的\(\sum\)裏，所以\(F(n)=\sum_{n|a}\sum_{n|b}\lfloor \frac Na\rfloor \lfloor \frac Mb \rfloor\)

所以答案為

\[f(1)=\sum_d\mu(d)\sum_{d|a}\sum_{d|b}\lfloor \frac Na\rfloor \lfloor \frac Mb \rfloor\]

發現後面的式子不好求，可以構造\(g(x)=\sum_{i=1}^x\lfloor \frac xi\rfloor\)

則最終的式子為

\[Ans=\sum_d\mu(d)g(\lfloor \frac nd\rfloor)g(\lfloor \frac md \rfloor)\]

分塊解決即可

至於如何求\(g(x)\)，可以發現其實\(g(n)=\sum_{i=1}^nd(i)\)，利用線性篩完約數函數\(d(i)\)後求個前綴和即可

Code

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

inline void read(int&x){
    char c11=getchar();x=0;while(!isdigit(c11))c11=getchar();
    while(isdigit(c11))x=x*10+c11-‘0‘,c11=getchar();
}

const int N=50101;
int pri[N],is[N];
int u[N],g[N],d[N];
int tp;ll sm;

int main(){
    g[1]=u[1]=1;
    for(int i=2;i<N;++i){
        if(!is[i])pri[++tp]=i,u[i]=-1,g[i]=2,d[i]=1;
        for(int j=1,k;j<=tp and (k=i*pri[j])<N;++j){
            is[k]=1;
            if(i%pri[j])u[k]-=u[i],g[k]=g[i]*2,d[k]=1;
            else {
                u[k]=0,d[k]=d[i]+1;
                g[k]=g[i]/(d[i]+1)*(d[i]+2);
                break;
            }
        }
        g[i]+=g[i-1],u[i]+=u[i-1];
    }
    int n,m,T;read(T);
    while(T--){
        read(n),read(m),sm=0ll;
        if(n>m)swap(n,m);
        for(int i=1,j;i<=n;i=j+1){
            j=min(n/(n/i),m/(m/i));
            sm+=1ll*(u[j]-u[i-1])*g[n/i]*g[m/i];
        }printf("%lld\n",sm);
    }return 0;
}

題解-SDOI2015 約數個數和