bzoj2301 莫比乌斯反演+分段优化

转载请注明出处，谢谢http://blog.csdn.net/bigtiao097?viewmode=contents

题意：

给定整数a,b,c,d,k，求$a\leq x\leq c,\ \ b\leq y\leq d(1\leq a,b,c,d,k\leq 5\times 10^4)$且gcd(x,y)=k的数对(x,y)有多少对.

思路：

这个题比hdu1695 复杂一点，hdu1695可以看一下这里
和hdu1695相比，由于没有了a=c=1这个条件，需要在莫比乌斯反演的基础上容斥一下
莫比乌斯反演：

$$g(n) = \sum_{n|d} f(d) \Longleftrightarrow f(n) = \sum _{n|d} \mu(\frac{d}{n})g(d)$$
设 $f(n)$代表 $\ \ \ gcd(x,y)=n\ \ \ 的（x，y）的对数$
设 $g(n) 代表 gcd(x,y)\%n =0 \ \ 的（x，y）的对数$
这样就有了 $$g(n) = \sum_{n|d} f(d)$$
然后利用反演公式，因为 $g(n)$非常的好求， $g(n) = \frac{c}{n}\times \frac{d}{n}$
根据上述的反演公式直接计算就可以了
关于容斥，假设solve(m,n)代表 $1\leq x \leq m,\ 1\leq y \leq m$的(x,y)的对数
则 $ans = solve(b,d) - solve(a,d) - solve(b,c) + solve(a,c)$

但是上述这些还不足以AC，会TLE
我们还需要用分段来优化一下枚举因子的倍数的过程

for(int i=1,nex;i<=a;i=nex+1)
{
     nex = min(a/(a/i),b/(b/i));
     res += 1LL *(mu[nex]-mu[i-1])*(a/i)*(b/i);
}

这里的mu数组是莫比乌斯函数的前缀和
需要注意的是，若n/i=t,则t是满足a*i<=n的a的最大值，则n/(n/i)就是满足商为n/i的i的最大值
上述代码中，$i$~$nex这一段中a/i的值相同，b/i也是$

---

具体代码如下：
Result：Accepted $\ \ \ $Memory:1728 kb $\ \ \ $Time :10772 ms

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int maxn=5e4+5;
bool vis[maxn];
int prime[maxn];
int mu[maxn];
int n,a,b,c,d,k,x;
ll ans;
void mobius()
{
    memset(vis,0,sizeof vis);
    mu[1] = 1;
    int tot = 0;
    for(int i=2;i<maxn;i++)
    {
        if( !vis[i] )
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j=0;j<tot;j++)
        {
            if(i * prime[j] >=maxn) break;
            vis[i * prime[j]] = true;
            if( i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
                mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i=2;i<maxn;i++)
        mu[i]+=mu[i-1];
}
ll solve(int a,int b)
{
    ll res = 0;
    if(a>b)swap(a,b);
    for(int i=1,nex;i<=a;i=nex+1)
    {
        nex = min(a/(a/i),b/(b/i));
        res += 1LL *(mu[nex]-mu[i-1])*(a/i)*(b/i);
    }
    return res;
}
int main()
{
    mobius();
    int T;
    cin>>T;
    while(T--)
    {
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        a--;
        c--;
        n = min(b,d);
        ans = solve(b/k,d/k)+solve(a/k,c/k)-solve(a/k,d/k)-solve(c/k,b/k);
        printf("%lld\n",ans);
    }
    return 0;
}