莫比乌斯反演合集

容斥原理

给一个长度为 $n$ 的未知数列 $a[...]$。已知任意 $a[i]\le max$ ，和一个数列 $h[i]$ ， $h[i]$ 表示数列 $a[...]$ 中值为 $i$ 的倍数的数的个数。

求：每个 $i\in [1,max]$，$ans[i]$。其中 $ans[i]$ 意为值恰好为 $i$ 的数的个数。

考虑容斥原理：

（其实这就是莫反性质的证明之一。）

$ans[1]=h(1)-h(2)-h(3)-h(5)-...-h(max)+h(2\ast 3)+h(2\ast 5)+...+h((max-1)\ast max)-...$

其中，$h(1)$也是全集的大小。

写得更规范一点：$ans(1)={\sum }_S(-1{)}^{|S|}h(p(S))$。其中 $S$ 是一个质数集合，集合中每个质数都不同。$p(S)$ 是这些质数的乘积。当 $S$ 为空集时（注意，这也是被允许的），$h(\varnothing )$ 定义为 $h(1)$，即全集。

观察这个式子，发现，$(-1{)}^{|S|}$ 其实就是 $\mu (p(S))$。

而当 $S$ 不满足定义，即 $S$ 中存在多个相同质数，此时 $\mu (p(S))=0$，不会被计入答案，那么我们就可以通过写成莫比乌斯函数，来扩展 $S$ 的定义。

于是原式可以转化为 $ans(1)={\sum }_{i=1}^{max}\mu (i)h(i)$。

同样的，我们考虑 $ans(x)$。我们肯定可以将所有 $a[...]$ 都除以 x ，然后重新计一套 $h[...]$ 。然后用同样的公式来求。但是显然这样太慢了。

其实可以直接类比：$ans(x)={\sum }_{i=1}^{\lfloor \frac{max}{x}\rfloor }\mu (i)h(xi)$。

如果不懂怎么类比，像上面那样用容斥原理写开就行。

直接考虑莫反的用法：

对于一些函数f(n)，如果很难直接求出它的值，而容易求出其 倍数和 或 约数和g(n)，那么可以通过莫比乌斯反演简化运算，求得f(n)的值

f(n) 即问题中的$ans[i]$，g(n) 即 $h[i]$ （倍数和）。

因此，可以直接得出：$h(n)={\sum }_{d:n|d}ans(d)$。

莫反： $g(n)={\sum }_{d:n|d}f(d)\leftrightarrow f(n)={\sum }_{d:n|d}g(d)\mu (\frac{d}{n})$

❗莫反从哪里来？

莫比乌斯函数的关键性质：${\sum }_{d|n}\mu (d)=\{\begin{array}{ll}1& n=1\\ 0& n\ne 1\end{array}$

即 ${\sum }_{d|n}\mu (d)=\epsilon (n)$, $\mu \ast 1=\epsilon$。（$\ast$为迪利克雷卷积）

因此，莫比乌斯函数是常数函数1的在迪利克雷卷积中逆元。一个函数与常数函数1进行迪利克雷卷积，得到的其实就是因数和或倍数和。那么我们要用因数和或倍数和反过来求这个函数，就直接卷上常数函数1的逆元，即莫比乌斯函数即可。

cf 990 G. GCD Counting

给一棵树，每个点一个权值，对于所有 $x\in (1,200000)$，问有多少条路径的 $gcd=x$。
考虑求一个更简单的问题：对于所有 $x\in (1,200000)$，问有多少条路经的 $gcd$ 是 $x$ 的倍数。
我们可以把原树分成200000张图。每张图都是原树拆掉特定点形成的。第 $i$ 张图中的点存在当且仅当这个点的值是 $i$ 的倍数，那么这棵树中的路径数就是对于这个 $i$ 我们要求的答案。路径数怎么算？对于每一个连通块，$C(n,2)$。但是我们要加上单个点。于是就成为了 $\sum \frac{n(n+1)}{2}$。（n为连通块的大小）
求出来更简单的问题 $h(x)$后，使用莫反即可得出原问题的答案。
$h(x)={\sum }_{x|d}ans(d)$
$ans(x)={\sum }_{x|d}\mu (\frac{d}{x})h(d)$
即
$ans(x)={\sum }_{i=1}^{\lfloor \frac{max}{x}\rfloor }\mu (i)h(xi)$

//
// Created by artist on 2021/9/18.
//


#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
#define all(x) x.begin(),x.end()
#define DB1(args...) do { cout << #args << " : "; dbg(args); } while (0)
#define endl '\n'
#define pii pair<int,int>

void dbg() { std::cout << "  #\n"; }

template<typename T, typename...Args>
void dbg(T a, Args...args) {
    std::cout << a << ' ';
    dbg(args...);
}

void io() { ios::sync_with_stdio(0), cin.tie(0), cout.tie(0); }

const int maxn = 200043;

int d[maxn];

void sieve(){
    vector<int> pr;
    for(int i=2;i<maxn;++i){
        if(!d[i]) d[i]=i,pr.pb(i);
        for(auto x:pr) {
            if(x>d[i]||x*1ll*i>=maxn) break;
            d[i*x] = x;
        }
    }
}

int mobius(int x){
    int lst = -1;
    int res = 1;
    while(x!=1){
        if(lst==d[x]) res=0;
        else res = -res;
        lst = d[x];
        x/=d[x];
    }
    return res;
}

int mb[maxn];

vector<int> need_bfs[maxn];
int used[maxn];
vector<int> g[maxn];
int a[maxn];
int cc = 0;
int q[maxn];
int hd,tl;
int good[maxn];

int bfs(int x) {
    hd = 0,tl = 0;
    q[tl++] = x;
    used[x] = cc;
    while(hd<tl) {
        int z=q[hd++];
        for(auto y:g[z]) {
            if(good[a[y]]==cc&&used[y]<cc) {
                used[y] = cc;
                q[tl++] = y;
            }
        }
    }
    return tl;
}

ll ans1[maxn];
ll ans2[maxn];

signed main() {
    //io();
    sieve();
    for(int i=1;i<maxn;++i) mb[i] = mobius(i);
    int n;cin>>n;
    for(int i=0;i<n;++i) cin>>a[i];
    for(int i=0;i<n-1;++i) {
        int x,y;cin>>x>>y;
        --x,--y;
        g[x].pb(y);
        g[y].pb(x);
    }
    for(int i=0;i<n;++i) need_bfs[a[i]].pb(i);
    for(int i=200000;i>=1;--i) {
        cc++;
        for(int j=i;j<=200000;j+=i)
            good[j] = cc;
        for(int j=i;j<=200000;j+=i)
            for(auto x:need_bfs[j]) {
                if(used[x]==cc) continue;
                int z = bfs(x);
                ans1[i] += z*1ll*(z+1)/2;
            }
        for(int j=i,k=1;j<=200000;j+=i,k++)
            ans2[i] += mb[k] * ans1[j];
    }
    for(int i=1;i<=200000;++i)
        if(ans2[i])
            cout<<i<<" "<<ans2[i]<<endl;
}

P1829 [国家集训队]Crash的数字表格 / JZPTAB

题意

求${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$。答案模20211009，n,m<=1e7。

推导

其实就是有两次d的枚举。单开一个函数提出来（简化问题）考虑最后再合并进去就好。

${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$

$={\sum }_{i=1}^n{\sum }_{j=1}^m\frac{ij}{\gcd (i,j)}$

$={\sum }_{d=1}{\sum }_{i=1}^n{\sum }_{j=1}^m\frac{ij}{d}[\gcd (i,j)=d]$

$={\sum }_{d=1}\frac{1}{d}{\sum }_{i=1}^n{\sum }_{j=1}^{m}ij[\gcd (i,j)=d]$

$={\sum }_{d=1}\frac{1}{d}{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }{d}^{2}ij[\gcd (i,j)=1]$

$={\sum }_{d=1}d{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }ij[\gcd (i,j)=1]$

设$g(n,m)={\sum }_{i=1}^n{\sum }_{j=1}^{m}ij[\gcd (i,j)=1]$

单考虑$g(n,m)$：

$g(n,m)={\sum }_{i=1}^n{\sum }_{j=1}^{m}ij{\sum }_{d|\gcd (i,j)}\mu (d)$

$={\sum }_{d=1}{d}^2\mu (d){\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }ij$

注意到${\sum }_{i=1}^n{\sum }_{j=1}^{m}ij=\frac{(n+1)n}{2}\times \frac{(m+1)m}{2}$，可以$O(1)$内求出。

因此$g(n,m)$可以分块求解。时间复杂度$O(\sqrt{N})$。

回到原式，发现也可以用分块求解。。。，于是答案就出来惹。

代码

//
// Created by acm on 2021/5/27.
//

#include<bits/stdc++.h>

using namespace std;
#define int long long int
#define mp(a, b) make_pair(a,b)
#define vi vector<int>
#define mii map<int,int>
#define mpi map<pair<int,int>,int>
#define vp vector<pair<int,int> >
#define pb(a) push_back(a)
#define fr(i, n) for(i=0;i<n;i++)
#define rep(i, a, n) for(i=a;i<n;i++)
#define F first
#define S second
#define endl "\n"
#define Endl "\n"
#define fast std::ios_base::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define mod 1000000007
#define dom 20101009
#define sl(a) (int)a.length()
#define sz(a) (int)a.size()
#define all(a) a.begin(),a.end()
#define pii pair<int,int>
#define mini 2000000000000000000
#define time_taken 1.0 * clock() / CLOCKS_PER_SEC

//mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
//const long double pi = acos(-1);
//mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count());
//primes for hashing 937, 1013
template<typename T, typename U>
static inline void amin(T &x, U y) {
    if (y < x) x = y;
}

template<typename T, typename U>
static inline void amax(T &x, U y) {
    if (x < y) x = y;
}

inline int read() {
    int x = 0, f = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9') {
        if (ch == '-') f = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9') {
        x = x * 10 + ch - '0';
        ch = getchar();
    }
    return x * f;
}

inline int qpow(int x, int y) {
    int res(1);
    while (y) {
        if (y & 1) res = 1ll * res * x % dom;
        x = 1ll * x * x % dom;
        y >>= 1;
    }
    return res;
}
const int maxn = 1e7+4;
int mu[maxn],sum[maxn],pri[maxn],cnt;
bool vis[maxn];

void getmu(int mx){
    mu[1]=1;
    for(int i=2;i<=mx;++i){
        if(!vis[i]) pri[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&pri[j]*i<=mx;++j){
            vis[i*pri[j]] = true;
            if(i%pri[j]==0){
                mu[i*pri[j]] = 0;
                break;
            }
            mu[pri[j]*i]=-mu[i];
        }
    }
    for(int i=1;i<=mx;++i) sum[i] = sum[i-1]+i*i%dom*(mu[i]+dom)%dom, sum[i]%=dom;
}

int solve3(int n,int m){
    return (n+1)*n/2%dom*((m+1)*m/2%dom)%dom;
}

int solve2(int n,int m){
    int ans = 0;
    for(int d=1,j;d<=min(n,m);d=j+1){
        j=min(n/(n/d),m/(m/d));
        ans += (sum[j]-sum[d-1]+dom)%dom*solve3(n/d,m/d)%dom;
        ans%=dom;
    }
    return ans;
}

int solve(int n,int m){
    int ans = 0;
    for(int d=1,j;d<=min(n,m);d=j+1){
        j=min(n/(n/d),m/(m/d));
        ans += (d+j)*(j-d+1)/2%dom*solve2(n/d,m/d)%dom;
        ans %= dom;
    }
    return ans;
}

signed main() {
    fast;
    int n,m;cin>>n>>m;
    getmu(max(n,m));
    cout<<solve(n,m)<<endl;
    return 0;
}