首先反演是学习orz vfk的，http://vfleaking.blog.uoj.ac/blog/87
补一个基础知识:怎么交换求和式
这里P(j)代表j满足的条件，并且[P()]表示当条件为真为1，否则值为0
关键的一个等式：
下面扯几道莫比乌斯反演的题目
bzoj2301 2045 1101 比较基础的莫比乌斯反演，首先拆成4个询问做，设f(i)表示在范围内i|gcd(a,b)的(a,b)对数，g(i)表示范围内gcd(a,b)=i的对数，则有
根据反演可得
我们知道最多只有
种取值，所以我们可以维护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.
bzoj2820 首先不难想到穷举质数，根据上面题的经验可得
多测这肯定炸掉，我们考虑优化，令T=pd则
不难发现，如果我们能求的前缀和，那么我们就变成了简单的bzoj2301
其实很简单，我们只要暴力穷举p更新求出所有h(T)即可，为什么，这里有两个定理
1.素数的个数近似n/lgn
2.
所以我们暴力更新的复杂度为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.
bzoj2154 首先，稍有基础的人都知道lcm(i,j)=ij/gcd(i,j)，所以
，我们令
,
则：
我们来进行反演
由开始推：
(公式恐惧症的不要跑……）
然后我们在处理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.
bzoj2693 2154的多测，我们考虑继续优化
下面我们只有求出后面那个（不妨叫作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.
bzoj3309 ……反正还是穷举最大公约数是吧，然后类似上题的设求、反演，我就不推了，最后得到
唯一的难点就是h(k)怎么快速求出，看起来最高指数幂是个奇怪的东西，但实际上我们设n=p1^a1*p2^a2……pn^an
如果任意ai=aj那么h(k)=(-1)^(n+1) 否则h(k)=0，这不(lan)难(de)证明