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=;
 int pri[mxn],mu[mxn],cnt=;
 bool vis[mxn];
 void init(){
     mu[]=;
     for(int i=;i<mxn;i++){
         if(!vis[i]){
             pri[++cnt]=i;
             mu[i]=-;
         }
         for(int j=;j<=cnt && pri[j]*i<mxn;j++){
             vis[pri[j]*i]=;
             if(i%pri[j]==){mu[pri[j]*i]=;break;}
             mu[pri[j]*i]=-mu[i];
         }
     }
     return;
 }
 LL calc(LL n){
     if(!n)return ;
     LL i,j,k,ed=floor(sqrt(n));
     LL res=,tmp=;
     for(k=;k<=ed;k++){
         if(mu[k]){
             tmp=;
             LL ED=n/(k*k);
             for(i=;i*i*i<=ED;i++){
                 for(j=i+;j*j*i<=ED;j++)
                     tmp+=(ED/(i*j)-j)*+;
                 tmp+=(ED/(i*i)-i)*;
                 tmp++;
             }
             res+=mu[k]*tmp;
         }
     }
     return (res+n)/;
 }
 LL a,b;
 int main(){
     init();
     scanf("%lld%lld",&a,&b);
     printf("%lld\n",calc(b)-calc(a-));
     return ;
 }