「算法笔记」莫比乌斯反演

莫比乌斯反演是数论中的一种小技巧,常用于解决一些数论函数的求和问题。通常做法是使用 μ μ 函数来将 [n=1] [ n = 1 ] 转换成 d|nμ(i) ∑ d | n μ ( i ) 的形式。

推荐博客:algocode 算法博客

μ μ 的定义:

  • 如果 n n 有至少一个平方因子,那么 μ(n)=0
  • 如果 n n 没有一个平方因子,且有 2k+1 个素因子,那么 μ(n)=1 μ ( n ) = − 1
  • 如果 n n 没有一个平方因子,且有 2k 个素因子,那么 μ(n)=1 μ ( n ) = 1

μ μ 的性质:

  • μ μ 是积性函数。
  • 根据这一性质,我们可以使用线性筛在 Θ(n) Θ ( n ) 的时间中筛出 i[1,n] i ∈ [ 1 , n ] 中的所有 μ(i) μ ( i )
void prework() {
    mu[1] = 1;
    for (int i = 2; i <= 50000; i++) {
        if (!mark[i]) {
            prime[++tot] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= tot && i * prime[j] <= 50000; j++) {
            mark[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
        mu[i] += mu[i - 1];
    }
}
  • d|nμ(i)=1 ∑ d | n μ ( i ) = 1
  • 证明见 algocode 算法博客。
  • 利用这一性质,我们就可以解决一些简单的问题了。

例题一:BZOJ 2301 Problem B

给定 a,b,c,d,e a , b , c , d , e ,求 bi=adj=c[gcd(i,j)=k] ∑ i = a b ∑ j = c d [ g c d ( i , j ) = k ] 的值。
1ab,cd,k50000 1 ≤ a ≤ b , c ≤ d , k ≤ 50000 50000 50000 组询问。

可以将询问拆成四个前缀询问。
问题变成了:求 ni=1nj=1[gcd(i,j)=k] ∑ i = 1 n ∑ j = 1 n [ g c d ( i , j ) = k ] 的值。

i,j i , j 都除以 k k ,上式等于 i=1nkj=1nk[gcd(i,j)==1]

利用 μ μ 函数的第二个性质,上式等于 nki=1nkj=1d|i,d|jμ(d) ∑ i = 1 ⌊ n k ⌋ ∑ j = 1 ⌊ n k ⌋ ∑ d | i , d | j μ ( d )

将枚举 d d 移到前面,并令 i=id,j=jd
则上式等于 nd=1μ(d)nkdi=1nkdj=11=nd=1μ(d)nkd2 ∑ d = 1 n μ ( d ) ∑ i ′ = 1 ⌊ n k d ⌋ ∑ j ′ = 1 ⌊ n k d ⌋ 1 = ∑ d = 1 n μ ( d ) ⌊ n k d ⌋ 2

ni ⌊ n i ⌋ 只有 Θ(n) Θ ( n ) 种取值。直接预处理 μ μ 函数的前缀和,然后数论分块即可。

#include <cstdio>
#include <iostream>
using namespace std;

const int maxn = 50005;
bool mark[maxn];
int tot, prime[maxn], mu[maxn];
int T, a, b, c, d, k;

void prework() {
    mu[1] = 1;
    for (int i = 2; i <= 50000; i++) {
        if (!mark[i]) {
            prime[++tot] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= tot && i * prime[j] <= 50000; j++) {
            mark[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
        mu[i] += mu[i - 1];
    }
}

int solve(int n, int m) {
    int res = 0;
    for (int i = 1, j; i <= min(n, m); i = j + 1) {
        j = min(n / (n / i), m / (m / i));
        res += (mu[j] - mu[i - 1]) * (n / i) * (m / i);
    }
    return res;
}

int main() {
    prework();
    scanf("%d", &T);
    while (T--) {
        scanf("%d %d %d %d %d", &a, &b, &c, &d, &k);
        printf("%d\n", solve(b / k, d / k) - solve(b / k, (c - 1) / k) - solve((a - 1) / k, d / k) + solve((a - 1) / k, (c - 1) / k));
    }
    return 0;
}

例题二:BZOJ 2154 Crash 的数字表格

求: ni=1mj=1lcm(i,j) ∑ i = 1 n ∑ j = 1 m l c m ( i , j ) n,m107 n , m ≤ 10 7

发现 lcm(i,j)=ijgcd(i,j) l c m ( i , j ) = i j g c d ( i , j )

考虑枚举 gcd(i,j)=d g c d ( i , j ) = d 。问题就变成了:

d=1min(n,m)di=1ndj=1md[gcd(i,j)=1]ij ∑ d = 1 m i n ( n , m ) d ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ [ g c d ( i , j ) = 1 ] ⋅ i j

f(n,m) f ( n , m ) 表示 ni=1mj=1[gcd(i,j)=1]ij ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = 1 ] ⋅ i j ,我们发现可以使用莫比乌斯反演得出原式等于:

d=1min(n,m)μ(d)i=1ndj=1ndij ∑ d = 1 m i n ( n , m ) μ ( d ) ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ i j

=d=1min(n,m)μ(d)nd(nd+1)2md(md+1)2 = ∑ d = 1 m i n ( n , m ) μ ( d ) ⌊ n d ⌋ ( ⌊ n d ⌋ + 1 ) 2 ⌊ m d ⌋ ( ⌊ m d ⌋ + 1 ) 2

这个式子可以使用数论分块来做。

然后,回到原问题。我们再使用一次数论分块,这样就可以在 Θ((n)2)=Θ(n) Θ ( ( n ) 2 ) = Θ ( n ) 的时间内求解了。

#include <cstdio>
#include <iostream>
using namespace std;
typedef long long ll;
const int maxn = 10000005;
const int mod = 20101009;
bool mark[maxn];
ll n, m, tot, prime[maxn], mu[maxn], sum[maxn];

void prework() {
    mu[1] = 1;
    for (int i = 2; i <= min(n, m); i++) {
        if (!mark[i]) {
            prime[++tot] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= tot && i * prime[j] <= min(n, m); j++) {
            mark[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (ll i = 1; i <= min(n, m); i++) {
        sum[i] = (sum[i - 1] + i * i % mod * (mu[i] + mod)) % mod;
    }
}

ll prod(ll x, ll y) {
    return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod;
}

ll func(ll x, ll y) {
    ll res = 0;
    for (ll i = 1, j; i <= min(x, y); i = j + 1) {
        j = min(x / (x / i), y / (y / i));
        res = (res + (sum[j] - sum[i - 1] + mod) * prod(x / i, y / i)) % mod;
    }
    return res;
}

int main() {
    scanf("%lld %lld", &n, &m);
    prework();
    ll ans = 0;
    for (ll i = 1, j; i <= min(n, m); i = j + 1) {
        j = min(n / (n / i), m / (m / i));
        ans = (ans + (j - i + 1) * (i + j) / 2 % mod * func(n / i, m / i)) % mod;
    }
    printf("%lld\n", ans);
    return 0;
}
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值