ZOJ3435 Ideal Puzzle Bobble

mobius反演的入门题。

首先，先抄一遍mobius反演：

已知 f(n) = sigma(d|n, g(d))

那么 g(n) = sigma(d|n, mu(d)*f(n/d))

mobius反演的另一种形式：

在某一范围内，已知 f(n) = sigma(n|d, g(d))

那么 g(n) = sigma(n|d, mu(d/n)*f(d))

本题的难度其实就入门来说比较大，推荐先做一下

POJ 3090, HDU1695, SPOJ VLATTICE

掌握本题之后，再去做更难一点的SPOJ PGCD

进入正题：

本题题意为求从(1,1,1)到(L,W,H)这个三维空间中的所有从(1,1,1)能看到的点，

换言之，求(0,0,0)到(L-1,W-1,H-1)的空间中所有(i,j,k)=1的点。

令 g(n) 为 (i,j,k) = n 的可见点的个数

f(n) 为 (i,j,k) | n 的可见点的个数

显然，我们可以得到  f(n) = sigma(n|d, g(d))

经过mobius反演后，可以推出 g(n) = sigma(n|d, mu(d/n)*f(d))

令 n = 1 ，则有g(1) = sigma(1|d, mu(d)*f(d)) 即本题要求的答案

函数f是可以直接计算的：

f(d) = (L/d) * (W/d) * (H/d)

所以g(1) = sigma(mu(d) * (L/d) * (W/d) * (H/d))

又因为本题是三维的情形，当需要统计的可见点在与(0,0,0)同一xy或xz或yz，将退化成二维情形，

故 ans = sigma(mu(d) * (L/d) * (W/d) * (H/d) + mu(d) * (L/d) * (W/d) + mu(d) * (L/d) * (H/d) + mu(d) * (H/d) * (W/d) )

= sigma(mu(d) * ( (L/d+1) * (W/d+1) * (H/d+1) - 1) )

再然后，本题是多case的卡时间，卡的特别狠！

除了预处理的时间之外，主程序中的时间复杂度即使是O(n)也过不了= =

于是只能继续优化了 = =

我们可以注意到 (n/d) 当d在某一连续的范围内时，值是不会变化的。

例如 (10/4) 和 (10/5) 就是一样的，(10/6) 到 (10/10) 都是一样的。

所以，我们可以预处理出前n项的mobius函数之和，然后就可以在O(1)的时间统计出这个区间内的mobius函数之和，再直接乘以这个区间任何一个 f 函数值就行了。

这样主程序的时间复杂度就会到O(sqrt(n))    （没错吧= =，好像是快一点点了= =）

至此，所有难点都已攻克。

怒贴代码：

#include <cstdio>#include <cstring>#include <algorithm>using namespace std;typedef long long ll;const int maxn = 1000001;int l,w,h;bool notprime[maxn];int prime[maxn],nprime;int mu[maxn],sum[maxn];ll ans;void getmu(){    nprime=0;    mu[1]=1;    for (int i=2;i<maxn;i++)    {        if (!notprime[i]) {prime[nprime++]=i; mu[i]=-1;}        for (int j=0;j<nprime&&i*prime[j]<maxn;j++)        {            notprime[i*prime[j]]=1;            if (i%prime[j]) mu[i*prime[j]]=-mu[i];            else {mu[i*prime[j]]=0; break;}        }    }    sum[0]=0;    for (int i=1;i<maxn;i++) sum[i]=sum[i-1]+mu[i];}int main(){    getmu();    while (scanf("%d%d%d",&l,&w,&h)==3)    {        l--; w--; h--;        ans=0;        for (int i=1,cur;i<=l||i<=w||i<=h;i=cur+1)        {            cur = 1000000;            if (i<=l) cur=min(cur, l/(l/i));            if (i<=w) cur=min(cur, w/(w/i));            if (i<=h) cur=min(cur, h/(h/i));            ll a=l/cur+1, b=w/cur+1, c=h/cur+1;            ans+=(sum[cur]-sum[i-1])*(a*b*c-1);        }        printf("%lld\n",ans);    }}