莫比乌斯反演入门

ACM算法日常

发布于 2020-06-22 10:37:23

1.1K0

发布于 2020-06-22 10:37:23

文章被收录于专栏：ACM算法日常

缘起

数论千万条, 反演第一条, 反演不会做, 队友两行泪... 所以我们来学一下莫比乌斯反演吧~

分析

《莫比乌斯函数入门》中我们已经学习了mobius函数. 现在要来进一步学习mobius反演.

首先思考一个有趣的问题

给定 N 和 M, 求满足 1<=x<=N, 1<=y<=M 且 gcd(x, y) 为质数的点对(x,y) 的数量.
数据范围 1<=N,M<=1e6

显然，暴力的 O(NM) 算法是绝对不能接受的.

此题的前置技能就是mobius反演.

何谓反演? 圆的反演

首先，我们先理解一下什么叫做反演? 反演的"演"的意思是演绎（deduction）的意思. 初中学习平面几何的时候就学过圆的反演公式. 如下图所示, 如果 , 其中 r 是圆的半径. 那么称为关于圆的反演点.

所以我们就知道什么叫做反演了——知道了的位置, 我们可以很轻松的定位的位置, 反之亦然.

反演是一种逻辑演绎的技巧, 如果一个问题（例如上面的）比较难处理，但是它的反演问题（例如上面的）却比较好处理, 则可以转换为其反演问题进行处理.

讲清楚了什么是反演，我们就来讲什么是mobius反演，为了讲清楚mobius反演就不得不提及dirichlet卷积.

Dirichlet 卷积

dirichlet卷积定义如下

显然，dirichlet卷积满足交换律和结合律

关于dirichlet卷积小例子如下

分别是不变映射和恒等映射, 那么

分别表示因子个数函数和因子和函数.

既然dirichlet卷积是一个运算，那么我们自然关心它的单位元，正如加法运算的 0 和乘法运算的1一样重要.

dirichlet卷积的单位元是

或者记做 , 其中 [] 我们称之为示性函数（数学中一般记做）.

这非常好证明. 只需要证明有即可

既然讲到了dirichlet卷积的单位元, 所以不得不提及dirichlet卷积的逆元. 事实上，我们一般关心mobius函数和欧拉函数的逆元. 而 nonvariant 函数就是的逆元. 也即

Proof: 按照定义, 只需要证明 , 我们有

而假设 n 有 k（k > 0）个不同的素因子. 则

证毕

至于欧拉函数，我们有

Proof: 和证明欧拉函数的计算公式类似. 令

我们证明 f 是一个积性函数. 事实上我们有

有了此积性公式，剩下的事情就很简单了. 令，则

其中用到了

所以我们就得到了

证毕.

一个值得注意的事实是, 上述欧拉函数的式子两边同时卷积上 , 我们有

展开写我们就得到了和这两个重要的积性函数之间的关系如下

好了, 介绍完毕dirichlet卷积，我们开始介绍mobius反演.

所谓mobius反演, 就是如下两个反演推导式子（回想一下本文开头对反演一词的解释）

站在反演的角度， g 是比较容易解决的问题, f 是比较难解决的问题.

上面两个式子的

相同点都是将比较难解决的 f 用比较容易解决的 g 表示了出来.
不同点是第一个式子是约数关系，第二个式子是倍数关系.

Proof:

第一个式子的证明

站在dirichlet卷积的角度, , 所以

其中，我们要注意到dirichlet卷积的交换律和结合律即可.

第二个式子的证明

令则考虑到的逆元是 nonvariant, 我们有

证毕.

至此，mobius反演讲解完毕. 我们看看如何使用mobius反演解决本文开头提出的问题.

首先，我们考虑如下定义

那么，显然的, 我们有

而且，幸运的是 , g 是一个非常好解决的问题

所以，我们利用mobius反演就可以计算 f 了

于是答案为

令

则答案进一步化简为

如果我们能 O(N) 预处理出 sum(d) 的话, 则本题的复杂度就是 O(N), 这是完全可以接受的.

事实上, 我们有如下预处理 sum 的伪代码

for (re i = 1; i <= tot; i++) // prime[1,...,tot] 是所有质数, m是考察上限
{
    for (re j = 1, t; (t = prime[i] * j) <= m; j++)
    {
        sum[t] += mu[j];
    }
}

于是就开心的解决了此问题了, 复杂度是 O(N) 的.

奉上完整的代码

// 此处省略一些头部文件...
const int maxn = 1e6+5;
int n, m, prime[maxn], isp[maxn], tot, mu[maxn], sum[maxn];

ilv mobius(int m)
{
 mu[1] = 1;
 for (re i = 2; i <= m; i++)
 {
  if (!isp[i])
  {
   prime[++tot] = i;
   mu[i] = -1;
  }
  for (re j = 1, t; j <= tot && (t = prime[j] * i) <= m; j++)
  {
   isp[t] = 1;
   if (i % prime[j])
   {
    mu[t] = -mu[i];
   }
   else
   {
    mu[t] = 0;
    break;
   }
  }
 }
}

signed main()
{
#ifdef LOCAL
 freopen("d:\\data.in", "r", stdin);
// freopen("d:\\my.out", "w", stdout);
#endif
 read(n), read(m);
 if (n > m)
 {
  swap(n, m);
 }
 mobius(m);
 for (re i = 1; i <= tot; i++)
 {
  for (re j = 1, t; (t = prime[i] * j) <= m; j++)
  {
   sum[t] += mu[j];
  }
 }
 int ans = 0;
 for (re i = 1; i <= n; i++)
 {
  ans += (n / i) * (m / i) * sum[i];
 }
 write(ans);
 flush();
 return 0;
}

应用上述程序可以快速计算出当 n = m = 1e6 时的答案为 274934018418，也就两千多亿——要是我有这么多钱就好了~

为了证明上面给出的mobius反演的正确性. 我们来看一道mobius反演的入门题 hdu 1695 gcd

给定 a, b, c, d, k, 求满足 a<=x<=b, c<=y<=d, (x,y)=k 的不同点对 (x, y) 的数量.
注意, (5, 7) 和 (7, 5) 被认为是相同点对

【输入】
多样例,第一行是样例数, 每个样例 有五个数.
0<a<=b<=1e5, 0<c<=d<=1e5, 0<=k<=1e5
输入保证 a = c = 1

【输出】
输出答案

【样例输入】
2
1 3 1 5 1
1 11014 1 14409 9

【样例输出】
Case 1: 9
Case 2: 736427

【hint】
对于第一个样例, 9对(x,y) 如下
(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 5), (3, 4), (3, 5).

【限制】
3s 32MB

首先，使用上面一样的推导，我们有

其中 f(1) 是集合 {} 的互质点对个数. 但是题目要求的是不同的互质点对的数量. 所以还要减去多算的点对. 而多算的点对恰好就是集合 {} 的互质点对个数的一半.

于是我们可以写出如下代码

//#include "stdafx.h"
//#define LOCAL
#pragma GCC optimize(2)
#pragma G++ optimize(2)
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <string>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
using namespace std;
#define int long long
#define re register int
#define ci const int
typedef pair<int, int> P;
#define FE(cur) for(re h = head[cur], to, len; ~h; h = g[h].nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt[cur].r - hjt[cur].l)
#define MID(cur) (hjt[cur].l + hjt[cur].r >> 1)
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
namespace fastio
{
 const int BUF = 1 << 21;
 char fr[BUF], fw[BUF], *pr1 = fr, *pr2 = fr;int pw;
 ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
 ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
 ilv pc(char c) { if (pw >= BUF) flush(); fw[pw++] = c; }
 ili read(int &x)
 {
  x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
  while(!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
  while(isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
  x *= f; return 1;
 }
 ili read(double &x) 
 {
  int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
  while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
  while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
  x = xx;
  if (c ^ '.') { x = f * xx; return 1; }
  c = gc();
  while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
  x *= f; return 1;
 }
 ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
 ilv writeln(int x) { write(x);pc(10); }
 ili read(char *x)
 {
  char c = gc(); if (!~c) return EOF;
  while(!isalpha(c) && !isdigit(c)) c = gc();
  while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
  *x = 0; return 1;
 }
 ili readln(char *x)
 {
  char c = gc(); if (!~c) return EOF;
  while(c == 10) c = gc();
  while(c >= 32 && c <= 126) *x++ = c, c = gc();
  *x = 0; return 1;
 }
 ilv write(char *x) { while(*x) pc(*x++); }
 ilv writeln(char *x) { write(x); pc(10); }
 ilv write(char c) { pc(c); }
 ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
const int maxn = 1e5+5;
int a, b, c, d, k, mu[maxn], isp[maxn], prime[maxn], tot;

ilv mobius(int m)
{
 mu[1] = 1;
 for (re i = 2; i <= m; i++)
 {
  if (!isp[i])
  {
   prime[++tot] = i;
   mu[i] = -1;
  }
  for (re j = 1, t; j <= tot && (t = i * prime[j]) <= m; j++)
  {
   isp[t] = 1;
   if (i % prime[j])
   {
    mu[t] = -mu[i];
   }
   else
   {
    mu[t] = 0;
    break;
   }
  }
 }
}

ili kk(int n, int m)
{
 int ans = 0;
 for (re i = 1; i <= n; i++)
 {
  ans += (n / i) * (m / i) * mu[i];
 }
 return ans;
}

signed main()
{
#ifdef LOCAL
 freopen("d:\\data.in", "r", stdin);
// freopen("d:\\my.out", "w", stdout);
#endif
 mobius(maxn - 1);
 int kase; scanf("%lld", &kase);
 for (re kse = 1; kse <= kase; kse++)
 {
  write("Case "), write(kse), write(": ");
  read(a), read(b), read(c), read(d), read(k);
  if (!k)
  {
   writeln("0");
   continue;
  }
  if (b > d) swap(b, d);
  b /= k, d /= k;
  writeln(kk(b, d) - kk(b, b) / 2);
 }
 flush();
 return 0;
}

ac情况如下

Accepted 93ms 5460kB

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-06-17，如有侵权请联系 cloudcommunity@tencent.com 删除

sum

本文分享自 ACM算法日常微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

sum

登录后参与评论

0 条评论

热度