【BZOJ2301】problem b，数论之莫比乌斯反演

Time:2016.05.27
Author:xiaoyimi
转载注明出处谢谢

---

传送门
思路：
$∑^{d}_{i=c}∑^{b}_{j=a}[gcd(i,j)=k]$

$=∑^{d}_{i=1}∑^{b}_{j=1}[gcd(i,j)=k]-∑^{c-1}_{i=1}∑^{b}_{j=1}[gcd(i,j)=k]-∑^{d}_{i=1}∑^{a-1}_{j=1}[gcd(i,j)=k]+∑^{c-1}_{i=1}∑^{a-1}_{j=1}[gcd(i,j)=k]$

可以观察到上面这四个式子都是类似的，我们只用处理一个就好
这里选用
$∑^{d}_{i=1}∑^{b}_{j=1}[gcd(i,j)=k]$

令$i=ki,j=kj$

$=∑^{⌊\frac d k⌋}_{i=1}∑^{⌊\frac b k⌋}_{j=1}[gcd(i,j)=1]$

$=∑^{⌊\frac d k⌋}_{i=1}∑^{⌊\frac b k⌋}_{j=1}∑_{p|gcd(i,j)}μ(p)$

设$b<=d,a<=c$

$=∑^{⌊\frac d k⌋}_{i=1}∑^{⌊\frac b k⌋}_{j=1}∑_{[p|i][p|j]}μ(p)$

$=∑_{p=1}^{⌊\frac b k⌋}⌊\frac b {pk}⌋⌊\frac d {pk}⌋μ(p)$

μ可以前缀和处理，然后$⌊\frac b {pk}⌋⌊\frac d {pk}⌋$的取值只有$\sqrt n$个，每次$O(n\sqrt n)$分块处理，总复杂度$O(n\sqrt n)$
注意：
1.实测不用开long long
2.推式子一定要细心，我刚开始就推错了好几次，尤其是分别处理4个式子，并且搞混了abcd……
3.好久不碰数论，欧拉筛都不会了……
代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#define M 50004
using namespace std;
int n,a,b,c,d,k;
int mu[M],prime[M];
bool vis[M];
int in()
{
    char ch=getchar();int t=0;
    while (!isdigit(ch)) ch=getchar();
    while (isdigit(ch)) t=(t<<3)+(t<<1)+ch-48,ch=getchar();
    return t;
}
int solve(int x,int y)
{
    if (x>y) swap(x,y);
    int ans=0,last;
    for (int i=1;i<=x;i=last+1)
    {
        last=min(x/(x/i),y/(y/i));
        ans+=(mu[last]-mu[i-1])*(x/i)*(y/i);
    }
    return ans;
}
main()
{
    mu[1]=1;
    for (int i=2;i<=50000;i++)
    {
        if (!vis[i])
            prime[++prime[0]]=i,
            mu[i]=-1;
        for (int j=1;j<=prime[0];j++)
            if (prime[j]*i>50000) break;
            else
            {
                vis[i*prime[j]]=1;
                if (i%prime[j])
                    mu[i*prime[j]]=-mu[i];
                else
                {
                    mu[i*prime[j]]=0;
                    break;
                }
            }
    }
    for (int i=2;i<=50000;i++) mu[i]+=mu[i-1];
    n=in();
    while (n--)
    {
        a=in()-1;b=in();c=in()-1;d=in();k=in();
        a/=k;b/=k;c/=k;d/=k;
        printf("%d\n",solve(b,d)-solve(b,c)-solve(d,a)+solve(a,c)); 
    }
}