【学习】容斥原理与莫比乌斯反演

容斥原理：

把判断一个东西的有无变成计算一个式子。

公式：

推论：

$$[n=0]=\sum_{i=0}^{n}(-1)^iC_{n}^{i}$$

证明：

用生成函数，有若干个个物品，考虑每一个选还是不选，假如说选择$n$个，可以得到：$(x-1)^n$。
其中当$x=1$的时候选择该物品，显然可以得到$[n==1]=(x-1)^n$。
用二项式定理展开左边，原式得证。

例题：

BZOJ 2440 [中山市选2011]完全平方数

题意

$T$组询问，每组输入一个$X$，问第$X$个不含平方因子的数是多少。（1不是完全平方数）
$1 \le X \le 10^9,T \le 50$

题解

我们可以先求出$1\to n$内不含平方因子数的个数，然后二分$n$的值。
现在的问题转化为求$1\to n$内不含平方因子数的个数，这个个数为$ans$。
可以对$n$以内有多少个质数构成因子进行容斥，$k$为质数个数的上限，$e$为质数的集合，$f(e)$为集合中元素的个数，$d$为元素$e$中的元素之积。

$$ans=\sum_{i=0}^{k}(-1)^i\sum_{e}^{f(e)=i}\lfloor \frac{n}{d^2} \rfloor$$
考虑到 $d$只有$\sqrt{n}$种取值，且 $(-1)^i=\mu(d)$，所以枚举 $d$。
$$ans=\sum^{\lfloor \sqrt{n}\rfloor}_{d=1}\mu(d)\lfloor \frac{n}{d^2} \rfloor$$

代码

#include <cstdio>
#include <cmath>
using namespace std;
const int maxn = 1000010, MX = 1e6;
typedef long long ll;
int mu[maxn], prime[maxn], pcnt;
bool is[maxn];
ll T, X;
void init(){
    is[0] = is[1] = 1, mu[1] = 1;
    for(register int i = 2; i <= MX; i ++){
        if(!is[i]) prime[++pcnt] = i, mu[i] = -1;
        for(register int j = 1; i*prime[j] <= MX; j ++){
            is[i*prime[j]] = 1;
            if(i%prime[j] == 0){mu[i*prime[j]] = 0; break;}
            else mu[i*prime[j]] = -mu[i];
        }
    }
}
inline ll calc(ll x){
    ll SZ = sqrt(x), ans = 0;
    for(register ll i = 1; i <= SZ; i ++){
        ans += mu[i] * (x/(i*i));
    }
    return ans;
}
int main(){
    scanf("%lld", &T);
    init();
    while(T --){
        scanf("%lld", &X);
        ll l = 0, r = 1e10, mid, ans = 0;
        while(l <= r){
            mid = ((l+r)>>1);
            if(calc(mid) < X) l = mid + 1;
            else ans = mid, r = mid - 1; 
        }
        printf("%lld\n", ans);
    }
    return 0;
}

BZOJ 2839 集合计数

题意

一个有$N$个元素的集合有$2^N$个不同子集（包含空集），现在要在这$2^N$个集合中取出若干集合（至少一个），使得它们的交集的元素个数为$K$，求取法的方案数，答案模$1e9+7$。
$1\le N\le 10^6,0\le K\le N$

题解

对除了$k$个以外的其他的交集个数进行容斥，设除了规定的$k$个交集，选出的其他集合的交集为$e$，设$f(e)$为$e$集合内元素的个数。
对$f(e)$进行容斥：

$$ans=C_{n}^{k}[f(e)=0]=C_{n}^{k}\sum_{i=0}^{n-k}(-1)^i \sum_{e}^{f(e)=i}C_{n-k}^{i}$$
现在的问题是计算大小为 $i$的交集集合$e$的可行的方案数，剩下的元素有 $n-k-i$个，每个考虑是否选择，可以得到 $2^{n-k-i}$个子集，在考虑每个子集是否选择，有 $2^{2^{n-k-i}}$种可能，因为要选出至少一个，所以答案是 $2^{2^{n-k-i}}-1$。
即：
$$ans=C_{n}^{k}\sum_{i=0}^{n-k}(-1)^iC_{n-k}^{i}(2^{2^{n-k-i}}-1)$$

代码

#include <cstdio>
#include <iostream>
using namespace std;
typedef long long ll;
const ll maxn = 1000010, p = 1e9+7;
ll jie[maxn], ni[maxn], ans, v, op = 1;
ll pow(ll x, ll y, ll pp){
    ll res = 1;
    for( ; y; y >>= 1, x = (x*x) % pp) if(y&1) res = (res * x) % pp;
    return res;
}
int n, k;
int main(){
    scanf("%d%d", &n, &k); jie[0] = jie[1] = ni[1] = ni[0] = 1;
    for(int i = 2; i <= n; i ++) jie[i] = (jie[i-1] * i) % p;
    for(int i = 2; i <= n; i ++) ni[i] = (((-(ll)(p/i)*ni[p%i])%p)+p) % p;
    for(int i = 2; i <= n; i ++) ni[i] = (ni[i]*ni[i-1]) % p;
    ans = (((jie[n] * ni[k]) % p) * ni[n-k]) % p;
    for(int i = 0; i <= n-k; i ++, op *= -1){
        ll pp = (((jie[n-k] * ni[n-k-i]) % p) * ni[i]) % p;
        ll t = pow(2, pow(2, n-k-i, p-1), p) - 1;
        v = (((v + op*(pp*t) % p) % p) + p) % p;
    }
    printf("%lld", (ans*v) % p);
    return 0;
}

莫比乌斯反演

POPOQQQ莫比乌斯反演课件
百度云链接 密码: rfqp

公式：

$$f(n)=\sum_{d|n}g(d) \Rightarrow g(n)=\sum_{d|n}\mu(d)f(\frac{n}{d})$$

推论：

$$[n=1]=\sum_{d|n}\mu(d)$$

例题：

BZOJ 2705 [SDOI2012]Longge的问题

题意

输入$n$，求

$$\sum_{i=1}^{n}gcd(i,n)$$
$0<N\le 2^{32}$

题解

和莫比乌斯反演没什么关系，主要是练习和式的计算。
令$g=gcd(i,n)$

$$ans=\sum_{g|n}g\sum_{i=1}^{n}[gcd(i,n)=g]=\sum_{g|n}g\sum_{i=1,g|i}^{n}[gcd(\frac{i}{g},\frac{n}{g})=1]$$
令 $i'=\frac{i}{g}$，所有的 $i$只能取$g$的倍数才可以，那 $i'$有意义的取值就是 $n$以内$g$的倍数的个数，所以：
$$ans=\sum_{g|n}g\sum_{i'=1}^{\frac{n}{g}}[gcd(i',\frac{n}{g})=1]$$
$$ans=\sum_{g|n}g\phi(\frac{n}{g})$$

代码

#include <cstdio>
#include <cmath>
typedef long long ll;
ll phi(ll x){
    ll up = sqrt(x);
    ll ans = x;
    for(int i = 2; i <= up; i ++){
        if(x%i == 0){
            ans = ans/i*(i-1);
            while(x%i == 0) x /= i;
        }
    }
    if(x != 1) ans = ans / x * (x - 1);
    return ans;
}
ll n, SQRT, res;
int main(){
    scanf("%lld", &n);
    SQRT = sqrt(n);
    for(int i = 1; i <= SQRT; i ++){
        if(n%i != 0) continue;
        res += i * phi(n/i);
        if(i*i != n) res += (n/i)*phi(i);
    }
    printf("%lld", res);
    return 0;
}

BZOJ 2301 [HAOI2011]Problem b

题意

对于给出的$n$个询问，每次求有多少个数对$(x,y)$，满足$a\le x\le b$，$c\le y\le d$，且$gcd(x,y) = k$，$gcd(x,y)$函数为$x$和$y$的最大公约数。
$1\le n\le 50000，1\le a\le b\le 50000，1\le c\le d\le 50000，1\le k\le 50000$

题解

用容斥把问题转化为求$x\in[1, n]$，$y\in[i,m]$时的数对$(x,y)$的个数。

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=k]=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[gcd(i,j)=1]=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|i,d|j}\mu(d)$$
然后考虑枚举 $d$，那么$i,j$一定是 $d$的倍数。
$$ans=\sum_{d}^{\min(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor)}\mu(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor$$
一种更清爽的写法，令 $\lfloor\frac{n}{k}\rfloor$为 $n'$，令 $\lfloor\frac{m}{k}\rfloor$为 $m'$。
$$ans=\sum_{d}^{\min(n',m')}\mu(d)\lfloor\frac{n'}{d}\rfloor\lfloor\frac{m'}{d}\rfloor$$
因为 $\lfloor\frac{n'}{d}\rfloor,\lfloor\frac{m'}{d}\rfloor$分别只有 $\sqrt{n'}$和 $\sqrt{m'}$种取值，所以我们可以枚举这些取值。
考虑 $\lfloor\frac{a}{i}\rfloor$在哪些范围内取值相同，已知 $i$的大小，可以$O(1)$计算出相同取值的末端。
公式： $last=a/(a/i)$。运行一下下面的代码就很清楚了。

#include <cstdio>
int main(){
    for(int i = 1, a = 100; i <= 100; i ++)
        printf(" -> %d\t->%d \n", i, a/i);
    return 0;
}

代码

#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
const int maxn = 50010, MX = 50005;
int T, a, b, c, d, k, ans;
int prime[maxn], mu[maxn], pcnt;
bool is[maxn];
void init(){
    is[1] = is[0] = mu[1] = 1;
    for(int i = 2; i <= MX; i ++){
        if(!is[i]) prime[++pcnt] = i, mu[i] = -1;
        for(int j = 1; i*prime[j] <= MX; j ++){
            is[i*prime[j]] = 1;
            if(i%prime[j] == 0){mu[i*prime[j]] = 0; break;}
            else mu[i*prime[j]] = -mu[i];
        }
    }
    for(int i = 1; i <= MX; i ++) mu[i] += mu[i-1];
}
int calc(int n, int m){
    n /= k, m /= k;
    if(n > m) swap(n, m);
    int ans = 0, last;
    for(int i = 1; i <= n; i = last+1){
        last = min(n/(n/i), m/(m/i));
        ans += (n/i)*(m/i)*(mu[last]-mu[i-1]);
    }
    return ans;
}
int main(){
    scanf("%d", &T);
    init();
    while(T --){
        scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        ans = calc(b, d) + calc(a-1, c-1) - calc(a-1, d) - calc(b, c-1);
        printf("%d\n", ans);
    } 
    return 0;
}

BZOJ 2820 YY的GCD

题意

$T$组数据，每组给定$N,M$，求$1\le x\le N$, $1\le y\le M$且$gcd(x, y)$为质数的$(x, y)$有多少对。
$T = 10^4$
$N, M \le 10^7$

题解

$$ans=\sum_{p\in prime}^{\min(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=p]$$
$$ans=\sum_{p\in prime}^{\min(n,m)}\sum_{i=1}^{\lfloor \frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{p}\rfloor}[gcd(i,j)=1]$$
$$ans=\sum_{p\in prime}^{\min(n,m)}\sum_{i=1}^{\lfloor \frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{p}\rfloor}\sum_{d|i,d|j}\mu(d)$$
$$ans=\sum_{p\in prime}^{\min(n,m)}\sum_{d}^{\min(\lfloor\frac{n}{p}\rfloor,\lfloor\frac{m}{p}\rfloor)}\mu(d)\lfloor\frac{n}{pd}\rfloor\lfloor\frac{m}{pd}\rfloor$$
令 $D=pd$，枚举 $D$的取值。
$$ans=\sum_{D}^{\min(n,m)}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{p|D,p\in prime}^{\min(n,m)}\mu(\frac{D}{p})$$
我们要想办法去维护后面这一项的前缀和， $n$以内的质数有$\frac{n}{\log n}$个，我们枚举每个质数的倍数，复杂度均摊 $\log n$，这样就可以在 $O(n)$的时间内预处理出后面这些项的前缀和了。

代码

#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
const int maxn = 10000010, MX = 1e7+5;
int T, N, M;
int prime[maxn], mu[maxn], pcnt;
long long s[maxn];
bool is[maxn];
void init(){
    is[1] = mu[1] = 1;
    for(int i = 2; i <= MX; i ++){
        if(!is[i]) prime[++pcnt] = i, mu[i] = -1;
        for(int j = 1; i*prime[j] <= MX; j ++){
            int x = i*prime[j]; is[x] = 1;
            if(i%prime[j] == 0){mu[x] = 0; break;}
            else mu[x] = -mu[i];
        }
    }
    for(int i = 1; i <= pcnt; i ++){
        for(int j = 1; j*prime[i] <= MX; j ++){
            s[j*prime[i]] += mu[j];
        } 
    }
    for(int i = 2; i <= MX; i ++) s[i] += s[i-1];
}
long long calc(int n, int m){
    long long ans = 0;
    int last;
    if(n > m) swap(n, m);
    for(long long i = 1; i <= n; i = last+1){
        last = min(n/(n/i), m/(m/i));
        ans += (n/i)*(m/i)*(s[last]-s[i-1]);
    }
    return ans;
}
int main(){
    init();
    scanf("%d", &T);
    while(T --){
        scanf("%d%d", &N, &M);
        printf("%lld\n", calc(N, M));
    }
    return 0;
}

BZOJ 2005 [Noi2010]能量采集

题意

$n$行，$m$列，一个坐标$(x, y)$来表示一个植物，植物生产能量。$(0, 0)$有一台能量汇集机器，在汇集的过程中有一定的能量损失。如果一棵植物与能量汇集机器连接而成的线段上有k棵植物，则能量的损失为$2k + 1$。现在要计算总的能量损失。
$1 \le n, m \le 10^5$

题解

不难发现，一个点$(x,y)$与$(0,0)$连接的直线经过的植物数量$k=gcd(x,y)$。

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)\times 2-1$$
$$ans=\sum_{g}^{\min(n,m)}(g\times2-1)\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=g]$$
$$ans=\sum_{g}^{\min(n,m)}(g\times2-1)\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}[gcd(i,j)=1]$$
进行莫比乌斯反演。
$$ans=\sum_{g}^{\min(n,m)}(g\times2-1)\sum_{i=1}^{\lfloor \frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor \frac{m}{g}\rfloor}\sum_{d|i,d|j}\mu(d)$$
考虑枚举 $d$。
$$ans=\sum_{g}^{\min(n,m)}(g\times2-1)\sum_{d}^{\min(\lfloor\frac{n}{g}\rfloor,\lfloor\frac{m}{g}\rfloor)}\mu(d)\lfloor\frac{n}{dg}\rfloor\lfloor\frac{m}{dg}\rfloor$$
令 $D=pd$，枚举 $D$的取值。
$$ans=\sum_{D}^{\min(n,m)}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{g|D}\mu(\frac{D}{g})(g\times2-1)$$
然后我们就可以 $n\log n$预处理后面一项的前缀和。前面可以用 $O(\sqrt{n}+\sqrt{m})$的分块做。

代码

#include <cstdio>
#include <iostream>
using namespace std;
const int maxn = 100010;
int prime[maxn], mu[maxn], pcnt, n, m, MX;
bool is[maxn];
long long f[maxn]; 
void init(){
    mu[1] = 1; is[1] = 1;
    for(int i = 2; i <= MX; i ++){
        if(is[i] == 0){prime[++pcnt] = i, mu[i] = -1;}
        for(int j = 1; i*prime[j] <= MX; j ++){
            int x = i*prime[j];
            is[x] = 1;
            if(i%prime[j] == 0){mu[x] = 0; break;}
            else mu[x] = -mu[i];
        }
    }
    for(long long i = 1; i <= MX; i ++){
        for(long long j = 1; i*j <= MX; j ++){
            f[i*j] += mu[j]*(i*2-1);
        }
    }
    for(int i = 1; i <= MX; i ++) f[i] += f[i-1];
}
long long calc(long long n, long long m){
    long long ans = 0;
    int last;
    if(n > m) swap(n, m);
    for(int i = 1; i <= n; i = last+1){
        last = min(n/(n/i), m/(m/i));
        ans += (n/i)*(m/i)*(f[last] - f[i-1]);
    }
    return ans;
}
int main(){
    scanf("%d%d", &n, &m);
    MX = max(n, m) + 5, init();
    printf("%lld", calc(n, m));
    return 0;
}

BZOJ 3529 [SDOI2014]数表

题意

有一张$n×m$的数表，其第$i$行第$j$列$(1\le i\le n，1\le j\le m)$的数值为能同时整除$i$和$j$的所有自然数之和。给定$a$，计算数表中不大于$a$的数之和。
$1\le N,m\le 10^5,1 \le Q \le 2×10^4$

题解

定义$f(i)$为$i$的因数之和。

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}f(gcd(i,j))[f(gcd(i,j))\le a]$$
令$g=gcd(i,j)$，枚举$g$。
$$ans=\sum_{g}^{\min(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}f(g)[f(g)\le a][gcd(i,j)=g]$$

$$ans=\sum_{g}^{\min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}f(g)[f(g)\le a][gcd(i,j)=1]$$
$$ans=\sum_{g}^{\min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}f(g)[f(g)\le a]\sum_{d|i,d|j}\mu(d)$$
上一题老套路，枚举 $d$。
$$ans=\sum_{g}^{\min(n,m)}f(g)[f(g)\le a]\sum_{d}^{\min(\lfloor\frac{n}{g}\rfloor,\lfloor\frac{m}{g}\rfloor)}\mu(d)\lfloor\frac{n}{dg}\rfloor\lfloor\frac{m}{dg}\rfloor$$
我们期望用前缀和去预处理一些量，然后就可以省掉一维的枚举，但现在的哪一项都不适合用前缀和维护，所以继续改变形式。
令 $D=dg$，枚举 $D$。
$$ans=\sum_{D}^{\min(n,m)}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{g|D}f(g)[f(g)\le a]\mu(\frac{D}{g})$$
这样就可以对前面的下取整分 $\sqrt{n}$块，每次查询后面那一块的前缀和，考虑到有 $a$的限制条件，所以把询问按$a$排序，找到 $f(g)$满足条件的 $g$，把所有$g$的倍数加上 $f(g)\times \mu(\frac{D}{g})$。
这就需要单点修改与查询区间的和，用树状数组即可。

代码

#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = 20010, maxm = 100010, MX = 1e5+5;
ll MOD;
struct node{
    int n, m, a, num;
    ll ans;
}e[maxn];
int Q, P = 1, pcnt;
int prime[maxm], mu[maxm], f[maxm], lowp[maxm], c[maxm];
ll s[maxm];
bool is[maxm];
inline bool cmp1(node a, node b){return a.a < b.a;}
inline bool cmp2(node a, node b){return a.num < b.num;}
inline bool cmp3(int a, int b){return f[a] < f[b];}
inline int lowbit(int x){return x&(-x);}
inline void update(int x, int V){
    for(register int i = x; i <= MX; i += lowbit(i))
        s[i] = (s[i]+V)&MOD;
}
inline ll ask(int x){
    ll res = 0;
    for(register int i = x; i >= 1; i -= lowbit(i))
        res = (res+s[i])&MOD;
    return res;
}
void init(){
    mu[1] = is[1] = is[0] = f[1] = 1;
    for(register int i = 2; i <= MX; i ++){
        if(!is[i]){
            prime[++pcnt] = i;
            lowp[i] = i;
            mu[i] = -1;
            f[i] = i+1;
        }
        for(register int j = 1; i*prime[j] <= MX; j ++){
            is[i*prime[j]] = 1;
            if(i%prime[j] == 0){
                lowp[i*prime[j]] = lowp[i] * prime[j];
                mu[i*prime[j]] = 0;
                if(i == lowp[i]) f[i*prime[j]] = f[i] + (i*prime[j]);
                else f[i*prime[j]] = f[prime[j]*lowp[i]] * f[i/lowp[i]];
                break;
            }else{
                lowp[i*prime[j]] = prime[j];
                mu[i*prime[j]] = -mu[i];
                f[i*prime[j]] = f[i] * f[prime[j]];
            }
        }
    }
} 
inline ll calc(int n, int m){
    ll res = 0;
    int last;
    if(n > m) swap(n, m);
    for(register int i = 1; i <= n; i = last+1){
        last = min(n/(n/i), m/(m/i));
        res = (res + (n/i)*(m/i) * (ask(last) - ask(i-1)) ) & MOD;
    }
    return res;
}
void renew(int n, int m, int a){
    if(n > m) swap(n, m);
    for(register int &i = P; f[c[i]] <= a; i ++)
        for(register int j = c[i]; j <= MX; j += c[i])
            update(j, f[c[i]]*mu[j/c[i]]);
}
inline int gt(){  
    char _ch;  
    int _num = 0, _op = 1, _ok = 0;  
    while(1){  
        _ch = getchar();  
        if(_ch == '-') _op *= -1;  
        else if(_ch >= '0' && _ch <= '9') _num = _num*10 + _ch - '0', _ok = 1;  
        else if(_ok) return _op * _num;  
    }  
}  
int main(){
    MOD = ((ll)1<<31), MOD -= 1;init();
    for(register int i = 1; i <= MX; i ++) c[i] = i;
    sort(c+1, c+1+MX, cmp3);scanf("%d", &Q);
    for(register int i = 1; i <= Q; i ++)
        e[i].n = gt(), e[i].m = gt(), e[i].a = gt(), e[i].num = i;
    sort(e+1, e+1+Q, cmp1);
    for(register int i = 1; i <= Q; i ++){
        renew(e[i].n, e[i].m, e[i].a);
        e[i].ans = calc(e[i].n, e[i].m);
    }sort(e+1, e+1+Q, cmp2);
    for(register int i = 1; i <= Q; i ++)
        printf("%lld\n", e[i].ans&MOD);
    return 0;
} 

BZOJ 2693 jzptab

题意

$T$组数据，每组两个正整数$N,M$，求

$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$$ 对 $10^8+9$取模的值。
$T\le 10^4, N,M\le 10^7$

题解

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}$$
$$ans=\sum_{g}^{\min(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=g]\frac{ij}{g}$$
提取出 $g$,让$i,j$只枚举 $g$的倍数，为了让$i\times j$的值不变，所以乘上 $g^2$。
$$ans=\sum_{g}^{\min(n,m)}g\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}[gcd(i,j)=1]ij$$
$$ans=\sum_{g}^{\min(n,m)}g\sum_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{g}\rfloor}ij\sum_{d|i,d|j}\mu(d)$$
和上面类似，我们枚举 $d$的值，这样原来的$i,j$就相当于枚举了 $d$的倍数，所以说把$d$从式子中提取出来，就可以变成从 $1$开始枚举了。
然后用求和公式就好啦。
$$ans=\frac{1}{4}\sum_{g}^{\min(n,m)}g\sum_{d}^{\min(\lfloor\frac{n}{g}\rfloor,\lfloor\frac{m}{g}\rfloor)}d^2\mu(d)(1+\lfloor\frac{n}{dg}\rfloor)\times(1+\lfloor\frac{m}{dg}\rfloor)\times \lfloor\frac{n}{dg}\rfloor\times \lfloor\frac{m}{dg}\rfloor$$
枚举一个 $D=dg$。
$$ans=\sum_{D}^{\min(n,m)}D\frac{1}{2}\lfloor\frac{n}{D}\rfloor(1+\lfloor\frac{n}{D}\rfloor)\times\frac{1}{2}\lfloor\frac{m}{D}\rfloor(1+\lfloor\frac{m}{D}\rfloor)\sum_{d|D}\mu(d)d$$
令 $$f(x)=\sum_{d|D}\mu(d)d　　g(x)=\mu(d)d$$
因为 $g(a)g(b)=ab\times\mu(a)\mu(b)=ab\times\mu(ab)=g(ab)$，所以 $g$是积性函数，根据狄利克雷卷积可知$f=g\cdot1$，所以 $f$也是积性函数，然后用线筛预处理前缀和。

代码

#include <iostream>
#include <cstdio>
using namespace std;
const long long maxn = 10000010, mod = 1e8+9;
long long prime[maxn], MX, f[maxn], lowp[maxn], pcnt, T, n, m;
bool is[maxn];
void init(){
    is[1] = 1, f[1] = 1;
    for(int i = 2; i <= MX; i ++){
        if(is[i] == 0){
            prime[++pcnt] = i;
            lowp[i] = i;
            f[i] = (1-i)%mod;
        }
        for(int j = 1; prime[j]*i <= MX; j ++){
            int x = i*prime[j];
            is[x] = 1;
            if(i%prime[j] == 0){
                lowp[x] = prime[j] * lowp[i];
                if(lowp[i] == i) f[x] = f[i];
                else f[x] = (f[x/lowp[x]] * f[lowp[x]])%mod;
                break;
            }else{
                lowp[x] = prime[j];
                f[x] = (f[i]*f[prime[j]])%mod;
            }
        }
    }
    for(int i = 2; i <= MX; i ++) f[i] = (f[i]*i%mod+f[i-1])%mod;
}
int main(){
    MX = 10000001; init();
    scanf("%d", &T);
    while(T --){
        scanf("%d%d", &n, &m);
        if(n > m) swap(n, m);
        long long last, ans = 0;
        for(long long i = 1; i <= n; i = last+1){
            last = min(n/(n/i), m/(m/i));
            ans = (ans + (((n/i)*(n/i+1)/2%mod)*((m/i)*(m/i+1)/2%mod)%mod)*(f[last] - f[i-1]))%mod;
        }
        printf("%lld\n", (ans+mod)%mod);
    }
    return 0;
}

BZOJ 3309 DZY Loves Math

题意

对于正整数$n$，定义 $f(n)$为 $n$所含质因子的最大幂指数。
例如$f(1960)=f(2^3 * 5^1 * 7^2)=3,f(10007)=1,f(1)=0$。
给定正整数 $a,b$，求

$$\sum_{i=1}^{a}\sum_{j=1}^{b}f(gcd(i,j)$$

题解

类似于上面的题目，很显然：

$$ans=\sum_{D}^{\min(n,m)}\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor\sum_{g|D}f(g)\mu(\frac{D}{g})$$
令 $$g(x)=\sum_{g|x}f(g)\mu(\frac{x}{g})$$
然后考虑如何求出 $g$函数的前缀和。
因为$\frac{x}{g}$再质因数的幂大于 $2$时$\mu(\frac{x}{g})$的值为0，对答案没有贡献，所以有贡献的 $f(g)$对应的函数值，只会等于 $f(x)$或 $f(x)-1$。即从每一个质因子中取出一个质数或者不取。
设 $x$有$n$个质因子,他们最大的幂 $a$有$m$个，分一下的情况进行讨论。

• 当$f(g)=f(x)$时。

  • 当存在$x$中有不同的质因数的时候，方案数为： $$ans=\sum_{i=0}^{m-1}(-1)^iC_{m}^{i}\sum_{j=0}^{n-m}(-1)^jC_{n-m}^{j}=0$$ 因为后式的值为$0$，所以这种情况下值恒为$0$。
  • 当$x$中所有的质因数全部相同时，若全部选择，那么方案数为： $$ans=\sum_{i=0}^{m}(-1)^iC_{m}^{i}=0$$
    既然最大的幂不可以都被取走，所以上我们现在要减去算错的方案，即： $$ans=-(-1)^m=(-1)^{m+1}$$

  这种情况下对答案的贡献为：$f(x)*ans=a\times(-1)^{m+1}$。

• 当$f(g)=f(x)-1$时。
  
  • 当有不同质因数的时候，方案数依旧为$0$。
  • 当所有质因数相同时，全部取走，对答案的贡献为$(a-1)\times(-1)^{m}$

综上所述，最终的答案为$-a\times(-1)^{m}+a\times(-1)^{m}-(-1)^{m}=(-1)^{m+1}$。
即我们求得：$g(x)=(-1)^{m+1}$，其中只有当$x$的质因数幂全部相同是对答案有贡献。
下一步我们考虑如何用线性筛法求$g$的前缀和，我们需要知道一个数除了最小质因数外，其他数的幂是否都相等，如果相等还需要知道等于多少，然后就是最小质因子的幂是多少，还有质因子的种类数，这样就可以筛出的值了。

代码

#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
const int maxn = 10000010;
typedef long long ll;
ll  g[maxn];
int l[10010], r[10010], prime[maxn], num[maxn], miy[maxn], lowp[maxn], mip[maxn];
int pcnt, n, m, MX, T;
bool is[maxn], ok[maxn];
void init(){
    is[1] = 1, g[1] = 0;
    for(int i = 2; i <= MX; i ++){
        if(is[i] == 0){
            lowp[i] = prime[++pcnt] = i;
            ok[i] = true; 
            num[i] = miy[i] = mip[i] = g[i] = 1;
        }
        for(int j = 1; prime[j]*i <= MX; j ++){
            int x = i*prime[j];
            is[x] = 1;
            if(i%prime[j] == 0){
                lowp[x] = lowp[i]*prime[j];
                mip[x] = mip[i] + 1, num[x] = num[i];
                if(lowp[x] == x) ok[x] = 1, miy[x] = miy[i]+1;
                else if(ok[x/lowp[x]] && mip[x/lowp[x]] == miy[x/lowp[x]]) ok[x] = 1, miy[x] = miy[x/lowp[x]];
                if(ok[x] && mip[x] == miy[x]) g[x] = (((num[x]+1)&1) ? -1 : 1);
                break;
            }else{
                lowp[x] = prime[j];
                mip[x] = 1, num[x] = num[i] + 1;
                if(ok[i] && mip[i] == miy[i]) ok[x] = 1, miy[x] = mip[i];
                if(ok[x] && mip[x] == miy[x]) g[x] = (((num[x]+1)&1) ? -1 : 1);
            }
        }
    }
    for(int i = 2; i <= MX; i ++) g[i] += g[i-1];
}
ll calc(ll n, ll m){
    ll res = 0, last;
    if(n > m) swap(n, m);
    for(ll i = 1; i <= n; i = last+1){
        ll cn=n/i,cm=m/i;
        last = min(n/cn, m/cm);
        res += cn*cm*(g[last]-g[i-1]);
    }
    return res;
}
inline int gt(){  
    char _ch;  
    int _num = 0, _op = 1, _ok = 0;  
    while(1){  
        _ch = getchar();  
        if(_ch == '-') _op *= -1;  
        else if(_ch >= '0' && _ch <= '9') _num = _num*10 + _ch - '0', _ok = 1;  
        else if(_ok) return _op * _num;  
    }  
}  
int main(){
    T = gt();
    for(int i = 1; i <= T; i ++){
        l[i] = gt(), r[i] = gt();
        if(l[i] > MX) MX = l[i];
        if(r[i] > MX) MX = r[i];    
    }
    init();
    for(int i = 1; i <= T; i ++) printf("%lld\n", calc(l[i], r[i]));
    return 0;
}