【HDU6134】Battlestation Operational (莫比乌斯反演)

多校第八场没想到碰到了一直想学却没有学的莫比乌斯反演，比赛时不出意外的打出了gg，便下定决心要好好会一会高深的莫比乌斯反演(在菜鸡的眼里是这样),下面是我的一些学习资料：

http://www.cnblogs.com/chenyang920/p/4811995.html

http://blog.csdn.net/acdreamers/article/details/8542292

https://wenku.baidu.com/view/fbec9c63ba1aa8114431d9ac.html

看完这些之后不得不佩服这些数学家的伟大，这个算法对于解决组合问题，包含gcd的条件的问题实在太厉害了。

然后就要开始讲这道题目了，这场的多校赛题目都是巨长，我就大概说一下这道题的题意：

题目给出一个n*n的矩阵对于矩阵中每一个点(x, y) 如果x >= y && gcd(x, y) = 1, 那么这个点的值为 x / y 向上取整， 否则就为0，要求这个矩阵的各个点的值的总和。即如公式所见：

分析：

首先设

因为gcd(i, j)等于1才有效，因此i / j 一定无法整除（除了j = 1），且结果就为 i / j + 1,因此有如下转化，其中φ(i)为欧拉函数

然后就是重要发现：

然后运用莫比乌斯反演

可得：

继续转化：

其中τ(n)表示n的约数个数，是积性函数也可以线性求出

其实后来发现：

F(n) = F（n - 1） + d(n - 1) + 1， 其中d(n)表示n的约数的个数

这样可以递推出每一项的F(n),

然后利用莫比乌斯反演公式就可以求出每一项f(n),注意这里求每一项的f(n)可以枚举j, 然后j的倍数的答案可以很容易求出，因为调和级数这样复杂度为nlog(n), f(n)的前缀和就是答案

#include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <map>#include <set>#include <vector>#include <queue>#include <cmath>#include <bitset>using namespace std;const int M = 1e9 + 7;const int MAXN = 1000005;int g[MAXN], pre[MAXN];bool check[MAXN+10];int prime[MAXN+10];int mu[MAXN+10];void Moblus(){    memset(check,false,sizeof(check));    mu[1] = 1;    int tot = 0;    for(int i = 2; i <= 1000000; i++)    {        if( !check[i] )        {            prime[tot++] = i;            mu[i] = -1;        }        for(int j = 0; j < tot; j++)        {            if(i * prime[j] > MAXN) break;            check[i * prime[j]] = true;            if( i % prime[j] == 0)            {                mu[i * prime[j]] = 0;                break;            }            else            {                mu[i * prime[j]] = -mu[i];            }        }    }}int factor[100][2];int getFactors(int x){    int fatCnt = 0;    int tmp = x;    for(int i = 0; prime[i] <= tmp / prime[i];i++)    {        factor[fatCnt][1] = 0;        if(tmp%prime[i] == 0)        {            factor[fatCnt][0] = prime[i];            while(tmp%prime[i] == 0)            {                factor[fatCnt][1]++;                tmp /= prime[i];            }            fatCnt++;        }    }    if(tmp != 1)    {        factor[fatCnt][0] = tmp;        factor[fatCnt++][1] = 1;    }    return fatCnt;}int main(){    Moblus();    g[1] = 1;    int cnt = 1;    for(int i = 2; i <= 1000000; ++i)    {        g[i] = g[i - 1] + cnt + 1;        int m = getFactors(i);        cnt = 1;        for(int j = 0; j < m; ++j)             cnt *= (factor[j][1] + 1);    }    for(int j = 1; j <= 1000000; ++j)        for(int k = j, t = 1; k <= 1000000; k += j, t++)            pre[k] += mu[t] * g[j];    for(int i = 2; i <= 1000000; ++i)        pre[i] = (pre[i] + pre[i - 1]) % M;    int n;    while(scanf("%d", &n) != EOF)    {        printf("%d\n", pre[n]);    }    return 0;}

然后这里还附一份我觉得更简洁的代码：

#include<cstdio>#include<cstring>#include<vector>using namespace std;typedef long long int ll;const ll mod = 1e9 + 7;const int A = 1e6 + 10;int mu[A],pri[A],tot;ll d[A],cnt[A],sum[A];bool vis[A];void init(){    tot = 0;mu[1] = d[1] = 1;    for(int i=2 ;i<A ;i++){        if(!vis[i]){mu[i] = -1;pri[++tot] = i;d[i] = 2;cnt[i] = 1;}        for(int j=1 ;j<=tot && pri[j]*i<A ;j++){            vis[i*pri[j]] = 1;            if(i%pri[j] == 0){                d[i*pri[j]] = d[i]/(cnt[i]+1)*(cnt[i]+2);                cnt[i*pri[j]] = cnt[i] + 1;                mu[i*pri[j]] = 0;                break;            }            d[i*pri[j]] = d[i]<<1;            cnt[i*pri[j]] = 1;            mu[i*pri[j]] = -mu[i];        }    }    sum[1] = 1;    for(int i=2 ;i<A ;i++){        sum[i] = (sum[i-1] + d[i-1] + 1)%mod;    }    for(int i=1 ;i<A ;i++){        sum[i] = (sum[i] + sum[i-1])%mod;        mu[i] = (mu[i] + mu[i-1])%mod;    }}int main(){    init();    int n;    while(~scanf("%d",&n)){        ll ans = 0;        int last;        for(int i=1 ;i<=n ;i=last + 1){            last = n/(n/i);            ans = (ans + (mu[last]-mu[i-1])%mod*(sum[n/i])%mod)%mod;        }        ans = (ans%mod + mod)%mod;        printf("%I64d\n",ans);    }    return 0;}