关于莫比乌斯反演的总结

首先反演是学习orz vfk的，http://vfleaking.blog.uoj.ac/blog/87

补一个基础知识:怎么交换求和式 $\sum_{P(j)} \sum_{P(k)} a_{j,k}=\sum_j \sum_k a_{j,k}\cdot [P(j,k)]=\sum_{P(k)} \sum_{P(j)} a_{j,k}$

这里P(j)代表j满足的条件，并且[P()]表示当条件为真为1，否则值为0

关键的一个等式：

下面扯几道莫比乌斯反演的题目

bzoj2301 2045 1101 比较基础的莫比乌斯反演，首先拆成4个询问做，设f(i)表示在范围内i|gcd(a,b)的(a,b)对数，g(i)表示范围内gcd(a,b)=i的对数，则有

$f(n)=\sum_{n|d} g(d)$ $f(n)=[\frac{x}{n}][\frac{y}{n}]$ 根据反演可得 $g(n)=\sum_{d} \mu (d) f(nd)=\sum_{d} \mu (d)\cdot [\frac{x}{nd}]\cdot [\frac{y}{nd}]$

我们知道 $[\frac{x}{d}]$ 最多只有 $2\sqrt{x}$ 种取值，所以我们可以维护miu的前缀和来搞定

 1 var a,b,c,d,k,t,i,j:longint;
 2     v:array[0..50005] of boolean;
 3     p,u:array[0..50005] of longint;
 4 
 5 function min(a,b:longint):longint;
 6   begin
 7     if a>b then exit(b) else exit(a);
 8   end;
 9 
10 procedure swap(var a,b:longint);
11   var c:longint;
12   begin
13     c:=a;
14     a:=b;
15     b:=c;
16   end;
17 
18 function cal(x,y:longint):int64;
19   var i,j:longint;
20   begin
21     if x>y then swap(x,y);
22     x:=x div k; y:=y div k;
23     cal:=0;
24     i:=1;
25     while i<=x do
26     begin
27       j:=min(x div (x div i),y div (y div i));
28       cal:=cal+int64(u[j]-u[i-1])*int64(x div i)*int64(y div i);
29       i:=j+1;
30     end;
31   end;
32 
33 begin
34   u[1]:=1;
35   for i:=2 to 50000 do
36   begin
37     if not v[i] then
38     begin
39       inc(t);
40       p[t]:=i;
41       u[i]:=-1;
42     end;
43     for j:=1 to t do
44     begin
45       if i*p[j]>50000 then break;
46       v[i*p[j]]:=true;
47       if i mod p[j]=0 then
48       begin
49         u[i*p[j]]:=0;
50         break;
51       end
52       else u[i*p[j]]:=u[i]*u[p[j]];
53     end;
54   end;
55   for i:=2 to 50000 do
56     u[i]:=u[i]+u[i-1];
57   readln(t);
58   while t>0 do
59   begin
60     dec(t);
61     readln(a,b,c,d,k);
62     writeln(cal(b,d)-cal(a-1,d)-cal(b,c-1)+cal(a-1,c-1));
63   end;
64 end.

2301

bzoj2820 首先不难想到穷举质数，根据上面题的经验可得

$ans=\sum_{p}\sum_{d} \mu (d)\cdot [\frac{x}{pd}]\cdot [\frac{y}{pd}]$ 多测这肯定炸掉，我们考虑优化，令T=pd则

$ans=\sum_{p}\sum_{p|T} \mu (\frac{T}{p})\cdot [\frac{x}{T}]\cdot [\frac{y}{T}]=\sum_{T}[\frac{x}{T}]\cdot [\frac{y}{T}]\sum_{p|T} \mu (\frac{T}{p})$

不难发现，如果我们能求 $h(T)=\sum_{p|T} \mu (\frac{T}{p})$ 的前缀和，那么我们就变成了简单的bzoj2301

其实很简单，我们只要暴力穷举p更新求出所有h(T)即可，为什么，这里有两个定理

1.素数的个数近似n/lgn

2. $\sum_{i=1}^n \frac{1}{i} = O(lgn)$

所以我们暴力更新的复杂度为O(n)

 1 const max=10000000;
 2 var u,f,p:array[0..max+1] of longint;
 3     v:array[0..max+1] of boolean;
 4     i,j,n,m,t,test:longint;
 5 
 6 function min(a,b:longint):longint;
 7   begin
 8     if a>b then exit(b) else exit(a);
 9   end;
10 
11 procedure swap(var a,b:longint);
12   var c:longint;
13   begin
14     c:=a;
15     a:=b;
16     b:=c;
17   end;
18 
19 function cal(x,y:longint):int64;
20   var i,j:longint;
21   begin
22     cal:=0;
23     if x>y then swap(x,y);
24     i:=1;
25     while i<=x do
26     begin
27       j:=min(x div (x div i),y div (y div i));
28       cal:=cal+int64(f[j]-f[i-1])*int64(x div i)*int64(y div i);
29       i:=j+1;
30     end;
31   end;
32 
33 begin
34   u[1]:=1;
35   for i:=2 to max do
36   begin
37     if not v[i] then
38     begin
39       inc(t);
40       p[t]:=i;
41       u[i]:=-1;
42     end;
43     for j:=1 to t do
44     begin
45       if i*p[j]>max then break;
46       v[i*p[j]]:=true;
47       if i mod p[j]=0 then
48       begin
49         u[i*p[j]]:=0;
50         break;
51       end
52       else u[i*p[j]]:=u[i]*u[p[j]];
53     end;
54   end;
55   for i:=1 to t do
56     for j:=1 to max div p[i] do
57       inc(f[p[i]*j],u[j]);
58   for i:=1 to max do
59     f[i]:=f[i]+f[i-1];
60   readln(test);
61   while test>0 do
62   begin
63     dec(test);
64     readln(n,m);
65     writeln(cal(n,m));
66   end;
67 end.

2820

bzoj2154 首先，稍有基础的人都知道lcm(i,j)=ij/gcd(i,j)，所以

$ans=\sum_d^{min(x,y)} \sum_{i}^{x} \sum_{j}^{y} \frac{ij}{d}\cdot [gcd(i,j)=d]$ ，我们令 $f(x,y)=\sum_{i}^x \sum_{j}^y ij$ , $g(x,y)=\sum_i^{x} \sum_j^{y} ij\cdot [i\perp j]$

则： $f(x,y)=\sum_{k}^{min(x,y)} {{k}^{2}}\cdot g([\frac{x}{k}],[\frac{y}{k}])$ $f(x,y)=\frac{(x+1)x}{2}\cdot \frac{(y+1)y}{2}$ $ans=\sum_d^{min(x,y)} \frac{g([\frac{x}{d}],[\frac{y}{d}])\cdot d^{2}}{d}=\sum_d^{min(x,y)}g([\frac{x}{d}],[\frac{y}{d}])\cdot d$

我们来进行反演

由 $\sum_{d|n} \mu (d)=[n=1]$ 开始推：

$g(x,y)=\sum_{m} [m=1]\cdot g([\frac{x}{m}],[\frac{y}{m}])\cdot {m}^{2}=\sum_{m} \sum_{k|m}\mu (k)\cdot g([\frac{x}{m}],[\frac{y}{m}])\cdot {m}^{2}=\sum_{k}\mu (k) \sum_{k|m}\cdot g([\frac{x}{m}],[\frac{y}{m}])\cdot {m}^{2}=\sum_{k}\mu (k) \sum_{m}^{mk\leq min(x,y)}\cdot g([\frac{x}{mk}],[\frac{y}{mk}])\cdot {(mk)}^{2}=\sum_{k}\mu (k) \cdot {k}^{2}\sum_{m}^{[\frac{min(x,y)}{k}]}\cdot g([\frac{[\frac{x}{k}]}{m}],[\frac{[\frac{y}{k}]}{m}])\cdot {m}^{2}=\sum_{k}\mu (k) \cdot {k}^{2}\cdot f([\frac{x}{k}],[\frac{y}{k}])$ (公式恐惧症的不要跑……）

然后我们在处理ans的时候是O(sqrt(x)),处理g(x,y)又是一个O(sqrt(x))，所以总的复杂度为O(x)

 1 const mo=20101009;
 2 var f:array[0..10000010] of int64;
 3     p,u:array[0..10000010] of longint;
 4     v:array[0..10000010] of boolean;
 5     i,j,n,m,t:longint;
 6     ans:int64;
 7 
 8 procedure swap(var a,b:longint);
 9   var c:longint;
10   begin
11     c:=a;
12     a:=b;
13     b:=c;
14   end;
15 
16 function min(a,b:longint):longint;
17   begin
18     if a>b then exit(b) else exit(a);
19   end;
20 
21 function sum(x:int64):int64;
22   var y:longint;
23   begin
24     y:=x+1;
25     if x mod 2=0 then x:=x div 2 else y:=y div 2;
26     exit(y*x mod mo);
27   end;
28 
29 function get(x,y:int64):int64;
30   begin
31     if x mod 2=0 then x:=x div 2 else y:=y div 2;
32     exit(x*y mod mo);
33   end;
34 
35 function cal(n,m:longint):int64;
36   var i,j:longint;
37   begin
38     i:=1;
39     cal:=0;
40     while i<=n do
41     begin
42       j:=min(n div (n div i),m div (m div i));
43       cal:=(cal+sum(n div i)*sum(m div i) mod mo*(f[j]-f[i-1]) mod mo) mod mo;
44       i:=j+1;
45     end;
46   end;
47 
48 begin
49   readln(n,m);
50   if n>m then swap(n,m);
51   u[1]:=1;
52   for i:=2 to n do
53   begin
54     if not v[i] then
55     begin
56       inc(t);
57       p[t]:=i;
58       u[i]:=-1;
59     end;
60     for j:=1 to t do
61     begin
62       if i*p[j]>n then break;
63       v[i*p[j]]:=true;
64       if i mod p[j]=0 then
65       begin
66         u[i*p[j]]:=0;
67         break;
68       end
69       else u[i*p[j]]:=-u[i];
70     end;
71   end;
72   for i:=1 to n do
73     f[i]:=(f[i-1]+int64(u[i])*int64(i)*int64(i) mod mo) mod mo;
74   i:=1;
75   while i<=n do
76   begin
77     j:=min(n div (n div i),m div (m div i));
78     ans:=(ans+get(j-i+1,j+i)*cal(n div i,m div i)) mod mo;
79     i:=j+1;
80   end;
81   writeln((ans+mo) mod mo);
82 end.

2154

bzoj2693 2154的多测，我们考虑继续优化

$ans=\sum_d g([\frac{x}{d}],[\frac{y}{d}])\cdot d=\sum_d\sum_{k}\mu (k) \cdot {k}^{2}\cdot f([\frac{x}{kd}],[\frac{y}{kd}])\cdot d$

下面我们只有求出后面那个（不妨叫作h(T)）的前缀和，就又变成bzoj2301拉

其实很简单，后面明显是积性函数，我们只要线性筛的时候求一下即可

具体的，当i mod p[j]=0 时 h(T)=h(i)*p[j] (因为miu(i*p[j])=0)，否则h(T)=h(i)*h(p[j])

 1 const mo=100000009;
 2       max=10000000;
 3 var f,u:array[0..10000010] of int64;
 4     p:array[0..10000010] of longint;
 5     v:array[0..10000010] of boolean;
 6     i,j,n,m,t,test:longint;
 7     ans:int64;
 8 
 9 procedure swap(var a,b:longint);
10   var c:longint;
11   begin
12     c:=a;
13     a:=b;
14     b:=c;
15   end;
16 
17 function min(a,b:longint):longint;
18   begin
19     if a>b then exit(b) else exit(a);
20   end;
21 
22 function sum(x:int64):int64;
23   var y:longint;
24   begin
25     y:=x+1;
26     if x mod 2=0 then x:=x div 2 else y:=y div 2;
27     exit(y*x mod mo);
28   end;
29 
30 function cal(n,m:longint):int64;
31   var i,j:longint;
32   begin
33     i:=1;
34     cal:=0;
35     while i<=n do
36     begin
37       j:=min(n div (n div i),m div (m div i));
38       cal:=(cal+sum(n div i)*sum(m div i) mod mo*(f[j]-f[i-1]) mod mo) mod mo;
39       i:=j+1;
40     end;
41   end;
42 
43 begin
44   u[1]:=1;
45   for i:=2 to max do
46   begin
47     if not v[i] then
48     begin
49       inc(t);
50       p[t]:=i;
51       u[i]:=(i-int64(i)*int64(i)) mod mo;
52     end;
53     for j:=1 to t do
54     begin
55       if i*p[j]>max then break;
56       v[i*p[j]]:=true;
57       if i mod p[j]=0 then
58       begin
59         u[i*p[j]]:=int64(p[j])*u[i] mod mo;
60         break;
61       end
62       else u[i*p[j]]:=u[i]*u[p[j]] mod mo;
63     end;
64   end;
65   for i:=1 to max do
66     f[i]:=(f[i-1]+u[i]) mod mo;
67   readln(test);
68   while test>0 do
69   begin
70     dec(test);
71     readln(n,m);
72     if n>m then swap(n,m);
73     writeln((cal(n,m)+mo) mod mo);
74   end;
75 end.

2693

bzoj3309 ……反正还是穷举最大公约数是吧，然后类似上题的设求、反演，我就不推了，最后得到

$ans=\sum_k [\frac{x}{k}][\frac{y}{k}]\cdot h(k)\,,h(k)=\sum_{d|k}\mu (\frac{k}{d}) f(d)$

唯一的难点就是h(k)怎么快速求出，看起来最高指数幂是个奇怪的东西，但实际上我们设n=p1^a1*p2^a2……pn^an

如果任意ai=aj那么h(k)=(-1)^(n+1) 否则h(k)=0，这不(lan)难(de)证明

posted on 2015-05-03 15:41 acphile 阅读(877) 评论(0) 编辑收藏举报