BZOJ_P2671 Calc(数论+莫比乌斯反演)

BZOJ传送门

Time Limit: 10 Sec Memory Limit: 128 MB
Submit: 192 Solved: 100
[Submit][Status][Discuss]
Description
给出N，统计满足下面条件的数对(a,b)的个数：
1.1<=a< b<=N
2.a+b整除a*b
　
Input
一行一个数N

　
Output
一行一个数表示答案

Sample Input
15

Sample Output
4

HINT
数据规模和约定
Test N Test N
1 <=10 11 <=5*10^7
2 <=50 12 <=10^8
3 <=10^3 13 <=2*10^8
4 <=5*10^3 14 <=3*10^8
5 <=2*10^4 15 <=5*10^8
6 <=2*10^5 16 <=10^9
7 <=2*10^6 17 <=10^9
8 <=10^7 18 <=2^31-1
9 <=2*10^7 19 <=2^31-1
10 <=3*10^7 20 <=2^31-1

Source

Sol:
拆出来公因数比较好想，但是跳表不太好打还有枚举约数也不太好弄，总之……我不会
题解
另:ans+=n/b/(a+b)*w;
这句话写成ans+=w*n/b/(a+b);
会WrongAnswer

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;
#define N 50005
//#define min(a,b) (a>b?b:a)
typedef long long ll;
ll top, stk[N];ll ans,n;
int mu[N],pr[N],cnt;bool b[N];
void GetMu(){
    mu[1]=1;
    for(int i=2;i<N;i++){
        if(!b[i]) pr[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*pr[j]<N;j++){
            b[i*pr[j]]=1;if(!(i%pr[j])) break;
            mu[i*pr[j]]=-mu[i];
        }
    }
}
void Decompose(ll x,ll i=1){
    for(top=0;i*i<x;i++) if(!(x%i)) stk[++top]=i,stk[++top]=x/i;
    if(i*i==x) stk[++top]=i;sort(stk+1,stk+top+1);
}
int main(){
    GetMu();scanf("%lld",&n);
    for(ll b=1;b*(b+1)<=n;b++){
        Decompose(b);
        for(ll a=1,j,w=0;a<b&&b*(a+b)<=n;a=j+1,w=0){
            j=min(n/(n/b/(a+b))/b-b,b-1);
            for(ll i=1;stk[i]<=j;i++) w+=mu[stk[i]]*(j/stk[i]-(a-1)/stk[i]);
            ans+=n/b/(a+b)*w;
        }
    }
    printf("%lld\n",ans);
    return 0;
}