[BZOJ2671][莫比乌斯反演]Calc

题目大意：给出$N$，求有多少组 $(a,b)$ 使得 $a+b|ab$ 。
我们首先令 $d=gcd(a,b)$ ， 则有 $a=id，b=jd$
$a + b | a b$
$=> d ( i + j ) | i j ( d ^ 2 )$
$=> ( i + j ) | i j d$
又因为 $gcd ( i , j ) == 1$
所以 $gcd ( i , i + j ) = gcd ( j , i + j ) = 1$
$=> ( i + j ) | d$
不妨设 $d = ( i + j ) k$，所以 $a = i ( i + j ) k , b = j ( i + j ) k$
即我们的问题现在转化为求有多少对三元组 $( i , j , k )$，
使得 $i ( i + j ) k <= N && j ( i + j ) k <= N && gcd ( i , j ) == 1$
若是枚举 $i < j$ 统计 k 的数量：$sigma( N / ( j ( i + j ) ) )$，时间复杂度 $O ( N )$ 。
对于 $N <= 2^31-1$ 我们是必然超时的，那么怎么优化呢？
我们会发现 $N / ( j ( i + j ) )$ 是整除，也就是说，
有一部分 $i , j$ 值的 $k$ 的数量是相同的，我们可以算出，
相同的 $k$ 的数量最多有 $2\sqrt N$ 种。
于是这里有一个神奇的做法：我们不妨枚举 $j ( j <= \sqrt N )$，
初始化 $i=1$ , 再令 $last = N / ( N / j / ( i + j ) ) / j - j$。
那么我们可以说 $[ i , last ]$ 这个区间中的 $N / ( j ( i + j ) )$ 是相同的。
至于上面说的神奇的方法，
可以理解为 本来 N 除以 $j ( i + j )$ 有余数可以使得在 $j$ 和 $N / ( j ( i + j ) )$ 不变的情况下，
i 还有增长的区间，除了两次后，$i$ 就变成了 $j$ 和 $N / ( j ( i + j ) )$ 不变时的的最大值。
这时下个区间就是$[ last + 1 , last' ]$
至于在当前区间内有多少个 $i$ 使得 $gcd(i,j)==1$ 就是莫比乌斯反演的基本内容了

具体细节可以看代码：

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdlib>
using namespace std;

#define MAXN 50000
#define MAXM
#define INF 0x3f3f3f3f
typedef long long int LL;

void Read(LL &x){
    x=0;char c=getchar();bool flag=0;
    while(c<'0'||'9'<c){if(c=='-')flag=1;c=getchar();}
    while('0'<=c&&c<='9'){x=x*10+c-'0';c=getchar();}
    if(flag)x=-x;
}

int u[MAXN+10];
bool isprime[MAXN+10];
int prime[MAXN+10],Cnt;

void init(int n){
    memset(isprime,1,sizeof(isprime));
    Cnt=0;

    n=sqrt(n);
    u[1]=1;
    for(int i=2;i<=n;++i){
        if(isprime[i]){
            prime[++Cnt]=i;
            u[i]=-1;
        }

        for(int j=1;j<=Cnt&&i*prime[j]<=n;++j){
            isprime[i*prime[j]]=0;
            if(i%prime[j])
                u[i*prime[j]]=-u[i];
            else{
                u[i*prime[j]]=0;
                break;
            }
        }
    }
}

int num[MAXN+10],top;
void Get_Divisor(LL n){
    LL i;
    top=0;
    for(i=1;i*i<n;i++)
        if(n%i==0)num[++top]=i,num[++top]=n/i;
    if(i*i==n)num[++top]=i;
    sort(num+1,num+top+1);
}

int main(){
    LL N;
    Read(N);
    init(N);

    LL i,j,k,last=0;
    LL cnt=0;
    LL ans=0;
    for(j=1;j*(j+1)<=N;++j){
        Get_Divisor(j);
        for(i=1;i<j&&j*(i+j)<=N;i=last+1){
            last=min(N/(N/j/(i+j))/j-j,j-1);
            cnt=0;
            for(k=1;num[k]<=last;++k)
                cnt+=u[num[k]]*(last/num[k]-(i-1)/num[k]);
            ans+=N/j/(i+j)*cnt;
        }
    }

    printf("%lld\n",ans);
}