莫比烏斯反演學習記錄(最菜的垃圾而淺薄基礎的總結)

鑑於容易忘,決定先把目前會的寫出來.....c++

莫比烏斯函數:
對於一個整數N,按照算數基本定理分解質因數爲N = p1^c1 * p2^c2 * p3^c3 * ... * pm^cm
              0 存在i∈[1, m],ci>1
μ(N) = { 1 m≡0(mod 2),任意i∈[1, m],ci = 1
            -1 m≡1(mod 2),任意i∈[1, m],ci = 1
通俗講:
①當N包含相等的質因子時,μ(N) = 0
②當N的全部質因子各不相等時,若N有偶數個質因子,μ(N) = 1
③當N有奇數個質因子時,μ(N) = -1
固然百度的話顯然能找到更好的說明方式,可是我以爲這樣比較清楚(雖然不利於理解求莫比烏斯函數)git

莫比烏斯函數求法markdown

 1 memset(vis, 0, sizeof(vis));  2     miu[1] = 1, vis[1] = 1;  3     for(int i = 2; i < v; ++i) {  4         if(!vis[i]) prime[++cnt] = i, miu[i] = -1;  5         for(int j = 1; j <= cnt; ++j) {  6             if(prime[j] > v / i) break;  7             vis[i * prime[j]] = 1;  8             if(i % prime[j] == 0) {//說明對於i * prime[j], 其中prime[j]的冪大於1 
 9                 miu[prime[j] * i] = 0; 10                 break; 11  } 12             miu[prime[j] * i] = -miu[i]; 13  } 14     }

 

莫比烏斯反演:
<1>莫比烏斯函數的一些性質:
·對於任意正整數n,Σμ(d) = [n=1] (d|n)
  證實:當n = 1時,顯然有函數g(n) = Σμ(d) = 1 (d|n)
     當n = p^k時(p爲質數)
      g(n) = Σμ(d) = μ(1)+μ(p)+μ(p^2)+...+μ(p^k)
             = 1 + (-1) + 0 + ... + 0 = 0
        當n = p1^c1 * p2^c2 * ... * pk^ck時
      g(n) = g(p1^c1)g(p2^c2)....g(pk^ck) = 0
·對於任意正整數n,Σμ(d)/d = φ(n)/n (d|n) //然而我並不會證實函數

<2>莫比烏斯反演
定理:F(n), f(n)爲兩個非負整數集合上的函數
   F(n) = Σf(d) (d|n) <==> f(n) = ΣΣμ(d)*f(e) (前一個Σ範圍爲d|n,後一個Σ範圍爲e|(n/d))
   證實:f(n)=Σμ(d)F(n/d) = Σμ(d)(d|n) * Σf(e)(e|(n/d)) = Σ(d|n)Σ(e|n/d)μ(d)*f(e)
                     交換求和順序有:f(n) = Σf(e)(e|n)Σμ(d)(d|n/e)
           根據莫比烏斯函數性質1,當且僅當n/e == 1,也就是n = e時,Σμ(d) = 1(d|n/e)
      則f(n) = f(e)*1 = f(n)
   證畢優化

(關於莫比烏斯反演,還有狄利克雷卷積證法,可是我不會.jpg)
莫比烏斯反演公式圖片搬運版(由於我根本不會markdown語法)ui

 

大佬眼中的入門題:
Bzoj2301spa

以cal(a, b, k)表示對於x<=a, y<=b,gcd(x, y)=k的x, y的對數
則咱們求的是
ΣΣ[gcd(i,j)==k](1<=i<=n, 1<=j<=m)
則根據容斥原理,答案顯然爲
cal(b, d, k) - cal(a - 1, d, k) - cal(b, c - 1, k) + cal(a - 1, b - 1, k)
/*
對於問題ΣΣ[gcd(i,j)==k](1<=i<=n, 1<=j<=m),能夠轉化爲
求對於1<=x<=[n/k],1<=y<=[m/k],gcd(x,y)互質的對數,也就是
ΣΣ[gcd(i,j)==1](1<=i<=[n/k], 1<=j<=[m/k])
實際上莫比烏斯反演式的本質就是Σμ(d)=[n==1] (d|n)
所以有Σμ(d)=[gcd(i,j)==1] (d|gcd(i,j))
所以也就是求
ΣΣΣμ(d) (1<=i<=[n/k], 1<=j<=[m/k], d|gcd(i,j))
.....
這些是另外一個常見推法的基本操做,可是我不會和式變換因此我不會這麼推....
*/
f(i)表示1<=x<=n, 1<=y<=m, 且gcd(x,y)=i的數對(x,y)的個數
F(i)表示1<=x<=n, 1<=y<=m, 且i|gcd(x,y)的數對(x,y)的個數
則F(i)=[n/i]*[m/i],
F[i]=Σf(d)(d|i)
→f(i)=Σμ(d/i)F(d) (i|d)
→f(i)=Σμ(d/i)[n/d][m/d] (i|d)
對於題目,咱們只須要求f(k),而後枚舉d(d爲k的倍數)
也就是枚舉k的每個倍數能夠在O(n)複雜度內回答一次詢問,顯然不能經過題目
考慮優化,容易注意到:咱們對於k的倍數的枚舉取決於[n/d]和[m/d]的取值
首先分析[n/d]的取值
1)當1<=d<=[sqrt(n)]時,最多有[sqrt(n)]種不一樣的取值
2)當[sqrt(n)]+1<=d<=n是,因爲[n/d]<=[sqrt(n)],所以有[sqrt(n)]種不一樣的取值
綜上所述+同理:
[n/d]有2[sqrt(n)]種不一樣的取值
[m/d]有2[sqrt(m)]種不一樣的取值
那麼[n/d][m/d]有多少種不一樣的取值呢
首先,不是4[sqrt(n)][sqrt(m)]個
考慮在一個數軸上,[n/d]將一段變成了2[sqrt(n)]個塊,每一個塊內取值相同,表示有2[sqrt(n)]中取值
同理考慮[m/d],而後手畫一個圖,能夠發現,對於[n/d]劃出的塊和[m/d]劃出的塊的間斷點是不一樣的
將[n/d]和[m/d]合併,實際上是間斷點的合併,即:塊的數量變爲了2[sqrt(n)]+2[sqrt(m)]
也就是說[n/d][m/d]一共有2[sqrt(n)]+2[sqrt(m)]個取值
對莫比烏斯函數維護前綴和,對於每一段就能夠直接回答了
令n <= m這樣對於每次詢問,複雜度爲O(sqrt(n))
總複雜度就是O(q*sqrt(n))3d

 1 #include<bits/stdc++.h>
 2 #define ll long long
 3 #define uint unsigned int
 4 #define ull unsigned long long
 5 using namespace std;  6 const int maxn = 50010;  7 int miu[maxn], sum_miu[maxn];  8 int T, a, b, c, d, k, cnt = 0;  9 int prime[maxn], vis[maxn]; 10 
11 inline int read() { 12     int x = 0, y = 1; 13     char ch = getchar(); 14     while(!isdigit(ch)) { 15         if(ch == '-') y = -1; 16         ch = getchar(); 17  } 18     while(isdigit(ch)) { 19         x = (x << 1) + (x << 3) + ch - '0'; 20         ch = getchar(); 21  } 22     return x * y; 23 } 24 
25 inline void pre_miu(int v) { 26     memset(vis, 0, sizeof(vis)); 27     miu[1] = 1, vis[1] = 1; 28     for(int i = 2; i < v; ++i) { 29         if(!vis[i]) prime[++cnt] = i, miu[i] = -1; 30         for(int j = 1; j <= cnt; ++j) { 31             if(prime[j] > v / i) break; 32             vis[i * prime[j]] = 1; 33             if(i % prime[j] == 0) {//說明對於i * prime[j], 其中prime[j]的冪大於1 
34                 miu[prime[j] * i] = 0; 35                 break; 36  } 37             miu[prime[j] * i] = -miu[i]; 38  } 39  } 40     for(int i = 2; i < v; ++i) miu[i] += miu[i - 1]; 41 } 42 
43 inline ll calc(int n, int m, int k) {//枚舉不一樣的值,last指向下一個間斷點,這個方法又被稱爲數論分塊 
44     n /= k, m /= k; 45     if(n > m) swap(n, m); 46     int last = 0; ll res = 0; 47     for(int i = 1; i <= n; i = last + 1) { 48         last = min(n / (n / i), m / (m / i)); 49         res += 1LL * (n / i) * (m / i) * (miu[last] - miu[i - 1]); 50  } 51     return res; 52 } 53 
54 int main() { 55     T = read(); 56  pre_miu(maxn); 57     while(T--) { 58         a = read(), b = read(), c = read(), d = read(), k = read(); 59         ll ans = calc(b, d, k) - calc(a - 1, d, k) - calc(b, c - 1, k) + calc(a - 1, c - 1, k); 60         printf("%lld\n", ans); 61  } 62     return 0; 63 }

後續的題目回頭慢慢補咕咕咕....code

Bzoj4804歐拉心算blog

 

對於AC此題,須要推式子(1)

關於此題的解法是如何想到的,我只能說我也是看課件(題解)的.....(2)

式子推的亂不要慌張,畢竟我不會markdown(3)

關於Σ的範圍和對應的字母,看我後面括號裏的順序肯定吧......(4)

對於式子ΣΣφ(gcd(i,j)) (1<=i<=n, 1<=j<=n),咱們令d = gcd(i, j)

可得ΣΣφ(d) (1<=i<=n, i<=j<=n),咱們進行和式變換,先枚舉d,可得:

  Σφ(d) ΣΣ[gcd(i, j) == d] (1<=d<=n , 1<=i<=n , 1<=j<=n)

==> Σφ(d) ΣΣ[gcd(i, j) == 1] (1<=d<=n , 1<=i<=[n/d] , 1<=j<=[n/d]) //注:此處以及如下的[x/y]形式的都表示下取整,(貌似數學內的取證符號就是下取整??(不清楚))

==> Σφ(d) Σμ(i)([n/(d*i)]^2)  (1<=d<=n , 1<=i<=[n/d]) //關於這一步如何推得?就是對後面兩個Σ進行反演

       令D = d*i, 和式變換得:

==> Σ([n/D]^2) Σφ(d)μ(D/d) (1<=D<=n, d|D)              //關於這一步如何推得?改變枚舉順序,將未知數只含D的項提出去,由於D=d*i,所以D∈[1, n],而後i用D/d表示,d的範圍就是d|D

這樣,接下來的問題轉變爲篩Σφ(d)μ(D/d) (d|D),觀察/打表/證實/猜測(或是看題解)可知這是一個積性函數

令h(d) = Σφ(d)μ(D/d) (d|D),考慮分類討論:設p爲一個質數

1) h(p) = φ(1)μ(p) + φ(p)μ(1) = -1 + p - 1 = p - 2

2) h(p^2) = φ(1)μ(p^2) + φ(p)μ(p) + φ(p^2)μ(1) = p^2 - 2p + 1 = (p - 1)^2

3) h(p^k) = φ(1)μ(p^k) + φ(p)μ(p^(k-1)) + ...... + φ(p^(k-1))μ(p) + φ(p^k)μ(1) = p^(k-1) (p-1) - p^(k-2) (p-1)

       = p^(k-2) (p-1)^2 = p^(k-2) h(p^2) = ph(p^k-1) //由此能夠獲得一個遞推式,設h(p^i), 則h(p^i) = h(p^(i-1))*p

4)gcd(a, p) = 1, 則根據積性函數, h(ap) = h(a)*h(p)

其他狀況根據積性函數推便可

 1 #include<bits/stdc++.h>
 2 #define ll long long
 3 #define uint unsigned int
 4 #define ull unsigned long long
 5 using namespace std;  6 const int maxn = 1e7+10;  7 int miu[maxn], prime[maxn];  8 bool vis[maxn];  9 int T, n, tot = 0; 10 int h[maxn]; 11 ll sum_h[maxn]; 12 
13 inline int read() { 14     int x = 0, y = 1; 15     char ch = getchar(); 16     while(!isdigit(ch)) { 17         if(ch == '-') y = -1; 18         ch = getchar(); 19  } 20     while(isdigit(ch)) { 21         x = (x << 1) + (x << 3) + ch - '0'; 22         ch = getchar(); 23  } 24     return x * y; 25 } 26 
27 inline void Mobius(int v) { 28     memset(vis, 0, sizeof(vis)); 29     h[1] = 1, vis[1] = 1; 30     for(int i = 2; i <= v; ++i) { 31         if(!vis[i]) { 32             vis[i] = 1; 33             prime[++tot] = i, h[i] = i - 2; 34  } 35         for(int j = 1; j <= tot; ++j) { 36             if(prime[j] > v / i) break; 37             vis[i * prime[j]] = 1; 38             if(i % prime[j] == 0) { 39                 if((i / prime[j]) % prime[j] == 0) 40                     h[i * prime[j]] = h[i] * prime[j]; 41                 else h[i * prime[j]] = h[i / prime[j]] * (prime[j] - 1) * (prime[j] - 1); 42                 //對於else 設i = ap, 則h(i*p) = h(ap^2) = h(a) * h(p^2),以及p^2 h(p^2) = 1 * (p-1)^2 
43                 break; 44  } 45             h[i * prime[j]] = h[i] * h[prime[j]];//無特殊狀況時,就是積性函數 
46  } 47  } 48     for(int i = 1; i <= v; ++i) sum_h[i] = sum_h[i - 1] + h[i]; 49 } 50 
51 int main() { 52  Mobius(maxn); 53     T = read(); 54     while(T--) { 55         n = read(); 56         ll ans = 0; 57         for(int i = 1, last; i <= n; i = last + 1) { 58             last = n / (n / i); 59             ans += 1LL * (n / i) * (n / i) * (sum_h[last] - sum_h[i - 1]); 60  } 61         printf("%lld\n", ans); 62  } 63     return 0; 64 }
相關文章
相關標籤/搜索