BZOJ2301 莫比乌斯反演

题意：a<=x<=b，c<=y<=d，求满足gcd(x,y)=k的数对(x,y)的数量         （(x,y)和(y,x)不算同一个）

比hdu1695多加了个下界，还有顺序不一样的也算上了。

因为G(x,y)本来就是顺序不一样的算不同方案，所以这题的公式就是：

Ans=G(b/k,d/k)-G((a-1)/k,d/k)-G(b/k,(c-1)/k)+G((a-1)/k,(c-1)/k)

但是本题数据很大，直接计算会TLE，

有一个优化：http://www.cnblogs.com/zhsl/p/3269288.html

//calc G(a,b)

LL _G(int a,int b)          //朴素算法
{
    int tx=min(a,b),ty=max(a,b);
    LL ans = 0;
    for(int i = 1; i <= tx; i++)
        ans += (LL)mu[i]*(tx/i)*(ty/i);
    return ans;
}

LL G(int n,int m)           //加分块优化
{
    LL ans = 0;
    if(n > m)   swap(n,m);
    for(int i = 1, la = 0; i <= n; i = la+1)
    {
        la = min(n/(n/i),m/(m/i));
        ans += (LL)(msum[la] - msum[i-1])*(n/i)*(m/i);      //事先预处理：msum[n]=SUM(mu[1..n])
    }
    return ans;
}

不知为什么自己测AC了结果bzoj上一直RuntimeError......要完蛋了orz，改成scanf和printf就过了

#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdio>
using namespace std;
#define LL long long
#define MMX 50100
LL mu[MMX],msum[MMX];
bool check[MMX];
int prime[MMX];
int a,b,c,d,k,T;

void Moblus()
{
    memset(check,false,sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i <= MMX; i++)
    {
        if( !check[i] )
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; j++)
        {
            if(i * prime[j] > MMX) break;
            check[i * prime[j]] = true;
            if( i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
            {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    msum[1]=mu[1];
    for (int i=2;i<=MMX;i++)
        msum[i]=msum[i-1]+mu[i];
}

//calc G(a,b)
LL _G(int a,int b)          //朴素算法
{
    int tx=min(a,b),ty=max(a,b);
    LL ans = 0;
    for(int i = 1; i <= tx; i++)
        ans += (LL)mu[i]*(tx/i)*(ty/i);
    return ans;
}

LL G(int n,int m)           //加分块优化
{
    LL ans = 0;
    if(n > m)   swap(n,m);
    for(int i = 1, la = 0; i <= n; i = la+1)
    {
        la = min(n/(n/i),m/(m/i));
        ans += (LL)(msum[la] - msum[i-1])*(n/i)*(m/i);      //事先预处理：msum[n]=SUM(mu[1..n])
    }
    return ans;
}

int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);

    cin>>T;
    Moblus();
    while (T--)
    {
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        //cout<<G(b/k,d/k)<<"  "<<G((a-1)/k,d/k)<<"  "<<G(b/k,(c-1)/k)<<"  "<<G((a-1)/k,(c-1)/k)<<endl;
        LL res=G(b/k,d/k)-G((a-1)/k,d/k)-G(b/k,(c-1)/k)+G((a-1)/k,(c-1)/k);
        printf("%lld
",res);
    }
    return 0;
}

实测：

朴素算法：   

1005

root

1002

Accepted

380K

47400MS

G++

1.6K

2014-11-15 15:37:03

优化算法：

1014

root

1002

Accepted

952K

6011MS

G++

1.81K

2014-11-15 17:17:59