[bzoj3994]约数个数和

前排提醒，本文是蒟蒻写的，极大可能性出错，请谨慎食用

好久不写数学题了，随便做一下吧。

题目大意：
给定N,M(<=50000),求：∑Ni=1∑Mj=1d(ij)

口胡：
首先我们有一个结论：∑Ni=1∑Mj=1d(ij)=∑Ni=1∑Mj=1⌊Ni⌋⌊Mj⌋[gcd(i,j)==1]
那么如何证明呢？
首先我们令F(N,M)=∑Ni=1∑Mj=1d(ij)
令G(N,M)=∑Ni=1∑Mj=1⌊Ni⌋⌊Mj⌋[gcd(i,j)==1]
做一个二维差分：

F(N,M)−F(N−1,M)−F(N,M−1)+F(N−1,M−1)=d(ij)

于是问题转化为了我们要证明：d(ij)=∑i|N∑j|M1∗[gcd(i,j)==1]
我们对于所有素数分类讨论，可以发现一个素数p，假设存在n=n′∗px,m=m′∗py那么对于式子左侧我们可以发现这个素数对于答案的贡献就是x+y+1，而式子右侧表示的其实就是1<=i<=N,1<=j<=m并且gcd(i,j)==1的有序数对的个数，那么我们考虑p，右边式子的贡献就是(px,1),(px−1,1).....(p,1),(1,1),(1,p),(1,p1).....(1,py)，也只有x+y+1种数对。

于是式子得到了证明：

∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)==1]

下面就是一个裸的反演啦~

∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)==1]=∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋∑d|gcd(i,j)μ(d)=∑d=1Nμ(d)∑i=1⌊Nd⌋∑j=1⌊Md⌋⌊Nid⌋⌊Mjd⌋=∑d=1Nμ(d)∑i=1⌊Nd⌋⌊Nid⌋∑j=1⌊Md⌋⌊Mjd⌋

我们令G(N)=∑Ni=1⌊Ni⌋
于是原式可变形为：∑Nd=1μ(d)G(⌊Nd⌋)G(⌊Md⌋)
那么G(n)函数就是群众喜闻乐见的d(N)的前缀和，线性筛约数个数即可。
线性筛复杂度：O(n)分块求和复杂度：O(n0.5)
code：

#include<bits/stdc++.h>using namespace std;const int maxn=5e4+10;typedef long long ll;ll prime[maxn],mu[maxn];bool isprime[maxn];ll c[maxn],d[maxn];int cnt;//c[i]表示i的最小因子的次数 inline void get(){    mu[1]=c[1]=d[1]=1;    for(int i=2;i<=maxn;i++){        if(!isprime[i]){            prime[++cnt]=i;            mu[i]=-1;            c[i]=1;            d[i]=2;        }        for(int j=1;j<=cnt&&i*prime[j]<maxn;j++){            isprime[i*prime[j]]=true;            if(i%prime[j]==0){                d[i*prime[j]]=d[i]*(c[i]+2)/(c[i]+1);                c[i*prime[j]]=c[i]+1;                break;            }            mu[i*prime[j]]=-mu[i];            d[i*prime[j]]=d[i]*d[prime[j]];            c[i*prime[j]]=1;        }    }    for(int i=1;i<maxn;i++)        mu[i]+=mu[i-1];    for(int j=1;j<maxn;j++)        d[j]+=d[j-1];}inline int read(){    int x=0;char ch=getchar();    while(!isdigit(ch)) ch=getchar();    while(isdigit(ch)) x=x*10+ch-48,ch=getchar();    return x;}int n,m;int main(int argc,const char * argv[]){    get();    int T=read();    while(T--){        n=read();m=read();        if(n>m) swap(n,m);        ll ans=0;        for(int i=1,nex;i<=n;i=nex+1){            nex=min(n/(n/i),m/(m/i));            ans+=(mu[nex]-mu[i-1])*d[n/i]*d[m/i];        }        printf("%lld\n",ans);    }    return 0;}