bzoj 3994 [SDOI2015]约数个数和

http://www.elijahqi.win/archives/3444
Description

设d(x)为x的约数个数，给定N、M，求

Input

输入文件包含多组测试数据。

第一行，一个整数T，表示测试数据的组数。
接下来的T行，每行两个整数N、M。

Output

T行，每行一个整数，表示你所求的答案。

Sample Input

2
7 4
5 6
Sample Output

110
121
HINT

1<=N, M<=50000

1<=T<=50000
Source

Round 1 感谢yts1999上传
考虑先化简原公式
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sigma(i*j)$
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{x|i}\sum\limits_{y|j}[gcd(x,y)==1]$
$\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}\sum\limits_{x|i}\sum\limits_{y|j}[gcd(x,y)==1]$
$\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}\sum\limits_{x|i}\sum\limits_{y|j}\sum\limits_{d|gcd(x,y)}\mu(d)$
$\sum\limits_{x=1}^{n}\sum\limits_{y=1}^{m}\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\sum\limits_{d|gcd(x,y)}\mu(d)$
$\sum\limits_{d=1}^{n}\mu(d)\sum\limits_{d|x}\sum\limits_{d|y}\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor$
$\sum\limits_{d=1}^{n}\mu(d)\sum\limits_{x=1}^{n}\sum\limits_{i=1}^{\frac{x}{d}}\lfloor\frac{n}{i\times d}\rfloor\sum\limits_{y=1}^{m}\sum\limits_{i=1}^{\frac{y}{d}}\lfloor\frac{m}{i\times d}\rfloor$
那么我们只要$n*\sqrt(n)$的预处理我后面那块内容即可
即$f[x]=\sum\limits_{i=1}^{x}\lfloor\frac{x}{i}\rfloor$
然后预处理莫比乌斯的前缀和即可在sqrt的时间内解决这题

#include<cstdio>
#include<cctype>
#include<algorithm>
#define ll long long
using namespace std;
inline char gc(){
    static char now[1<<16],*S,*T;
    if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(!isdigit(ch)) {if(ch=='-') f=-1;ch=gc();}
    while(isdigit(ch)) x=x*10+ch-'0',ch=gc();
    return x*f;
}
const int N=50050;
int f[N],prime[N],tot,mu[N],T,n,m;
bool not_prime[N];
inline void init(){
    for (int i=1;i<=5e4;++i){
        int last;
        for (int j=1;j<=i;j=last+1){
            last=i/(i/j);
            f[i]+=i/j*(last-j+1);
        }
    }mu[1]=1;
    for (int i=2;i<=5e4;++i){
        if (!not_prime[i]) prime[++tot]=i,mu[i]=-1;
        for (int j=1;prime[j]*i<=5e4;++j){
            not_prime[prime[j]*i]=1;
            if (i%prime[j]==0){
                mu[prime[j]*i]=0;break;
            }else mu[prime[j]*i]=-mu[i];
        }
    }
    for (int i=1;i<=5e4;++i) mu[i]=mu[i-1]+mu[i];
}
int main(){
//  freopen("bzoj3994.in","r",stdin);
    T=read();init();
    while(T--){
        n=read();m=read();int last;ll ans=0;
        for (int i=1;i<=min(n,m);i=last+1){
            last=min(n/(n/i),m/(m/i));
            ans+=(ll)(mu[last]-mu[i-1])*f[n/i]*f[m/i];
        }printf("%lld\n",ans);
    }
    return 0;
}