51nod1222 最小公倍数计数

题目来源： Project Euler

基准时间限制：6 秒 空间限制：131072 KB 分值: 640 

定义F(n)表示最小公倍数为n的二元组的数量。

即：如果存在两个数（二元组）X，Y(X <= Y)，它们的最小公倍数为N，则F(n)的计数加1。

例如：F(6) = 5，因为[2,3] [1,6] [2,6] [3,6] [6,6]的最小公倍数等于6。

 

给出一个区间[a,b]，求最小公倍数在这个区间的不同二元组的数量。

例如：a = 4，b = 6。符合条件的二元组包括：

[1,4] [2,4] [4,4] [1,5] [5,5] [2,3] [1,6] [2,6] [3,6] [6,6]，共10组不同的组合。

 

Input

输入数据包括2个数：a, b，中间用空格分隔（1 <= a <= b <= 10^11)。

Output

输出最小公倍数在这个区间的不同二元组的数量。

Input示例

4 6

Output示例

10

数学问题 莫比乌斯反演

请开始你的反演！
设：

$$ans(n)=sum_{i=1}^{n} sum_{j=1}^{n} [frac{i*j}{gcd(i,j)}<=n]$$
那么 $ans(b)-ans(a-1)$ 就是最终答案

尝试化简上面的式子：
$$sum_{i=1}^{n} sum_{j=1}^{n} [frac{i*j}{gcd(i,j)}<=n]$$
$$sum_{d=1}^{n} sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{n}{d}} [i*j<=frac{n}{d}] [gcd(i,j)==1]$$
$$sum_{d=1}^{n} sum_{k=1}^{frac{n}{d}} mu(k) sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{n}{d}} [i*k*j*k<=frac{n}{d}] $$
$$sum_{k=1}^{n} mu(k) sum_{d=1}^{frac{n}{k}} sum_{i=1}^{frac{n}{dk}} sum_{j=1}^{frac{n}{dk}} [i*j*d<=frac{n}{k^2}] $$

　　显然d和k值大到一定程度，最后面就是0了，所以我们可以缩小求和上界：

$$sum_{k=1}^{sqrt n} mu(k) sum_{d=1}^{frac{n}{k^2}} sum_{i=1}^{frac{n}{dk^2}} sum_{j=1}^{frac{n}{dk^2}} [i*j*d<=frac{n}{k^2}] $$

　　这个范围很友好，我们可以枚举$mu(k)$，求满足条件的i j d三元组数量。
　　需要求的三元组是无序的，为了不重不漏地计数，我们可以分别求出有序（单调上升）的三元组数量，对于其中三个数各不同的、有两个数相同的、三个数都相同的分别计数，然后乘以对应的组合数即可。

#include<iostream>
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<cmath>
#define LL long long
using namespace std;
const int mxn=1000010;
int pri[mxn],mu[mxn],cnt=0;
bool vis[mxn];
void init(){
    mu[1]=1;
    for(int i=2;i<mxn;i++){
        if(!vis[i]){
            pri[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt && pri[j]*i<mxn;j++){
            vis[pri[j]*i]=1;
            if(i%pri[j]==0){mu[pri[j]*i]=0;break;}
            mu[pri[j]*i]=-mu[i];
        }
    }
    return;
}
LL calc(LL n){
    if(!n)return 0;
    LL i,j,k,ed=floor(sqrt(n));
    LL res=0,tmp=0;
    for(k=1;k<=ed;k++){
        if(mu[k]){
            tmp=0;
            LL ED=n/(k*k);
            for(i=1;i*i*i<=ED;i++){
                for(j=i+1;j*j*i<=ED;j++)
                    tmp+=(ED/(i*j)-j)*6+3;
                tmp+=(ED/(i*i)-i)*3;
                tmp++;
            }
            res+=mu[k]*tmp;
        }
    }
    return (res+n)/2;
}
LL a,b;
int main(){
    init();
    scanf("%lld%lld",&a,&b);
    printf("%lld
",calc(b)-calc(a-1));
    return 0;
}