跳转至

Schreier–Sims 算法

引入

Schreier–Sims 算法 是計算羣論的一種算法,以數學家 Otto Schreier 和 Charles Sims 的名字命名。它可以在多項式時間(polynomial time)內找到有限置換羣的階數、查看給定排列是否包含在羣中以及其他許多任務。Schreier–Sims 算法最早由 Sims 在 1970 年基於 Schreier 的子羣引理引入。1991 年,Donald Knuth 基於此改進了時序。後來,又有了該算法更快的隨機版本。

註釋

羣論(Group theory)是數學的一個分支。在數學和抽象代數中,羣論研究稱為羣的代數結構。羣是由一組元素和一個可以應用於該集合的兩個元素的二元運算組成的系統。

背景

註釋

具體算法和時間複雜度分析請看之後章節,此節只是簡單的背景介紹。

Schreier–Sims 算法是一種計算置換羣的基強生成集(BSGS, base and strong generating set)的有效方法。特別是,SGS 確定了羣的順序並使查看給定排列是否包含在羣中變得更加容易。由於 SGS 對於計算羣論中的許多算法至關重要,因此計算機代數系統通常依賴 Schreier–Sims 算法進行羣的有效計算。

Schreier–Sims 的運行時間因實現而異。令 \(G\leq S_{n}\)\(t\) 個生成器給出。該算法可能的運行時間為:

  • \(O(n^{2}\log ^{3}|G|+tn\log |G|)\) 需要 \(O(n^{2}\log |G|+tn)\) 的內存;
  • \(O(n^{3}\log ^{3}|G|+tn^{2}\log |G|)\) 需要 \(O(n\log ^{2}|G|+tn)\) 的內存;

Schreier 向量的使用會對 Schreier–Sims 算法的實現性能產生重大影響。對於 Schreier–Sims 算法的 Monte Carlo 變體,預估複雜度如下:

  • \(O(n\log n\log ^{4}|G|+tn\log |G|)\) 需要 \(O(n\log |G|+tn)\) 的內存;

現代計算機代數系統(例如 GAP 和 Magma)通常使用優化過的蒙特卡羅算法。

定義

定義任意羣 \(G\) 是由一個子集 \(S\) 生成,\(G = \langle S \rangle\), 對於所有 \(i\),\(G\) 的每個元素都可以寫成 \(s_{1}, \cdots, s_{r}\)\(s_{i} \in S\)\(s_{i}^{-1} \in S\) 的乘積。

Schreier 樹

對於 \(S\) 根為 \(\alpha\)Schreier 樹\(α\) 軌道的如下表示:

  • Schreier 樹是一棵以 \(\alpha\) 為根,以 \(\alpha^{G}\) 的元素為頂點的樹,
  • 它的邊描述了從 \(\alpha\) 到每個頂點所需的 \(S\) 的元素,即樹中的每條邊 \({i, j}\),其中 比 \(j\) 更靠近根的 \(i\) 由生成器 \(s \in S\) 標記,將 \(i\) 移動到 \(j\)
  • Schreier 樹可以通過 廣度優先搜索深度優先搜索\(\alpha\) 開始用所有生成器 \(s \in S\) 嘗試到達新的節點 \(\alpha^{s}\) 來找到。因此,計算 Schreier 樹所需的時間以 \(O(rn)\) 為界。這意味着我們可以以一種有效的方式找到 \(|\alpha^{G}|\)

Schreier 引理

Schreier 引理説明了如何使用 \(\alpha\) 的 Schreier 樹來查找穩定器 \(G_{\alpha}\) 的生成器。

Schreier 引理:令 \(G = \langle S \rangle\)。那麼由一組 Schreier 生成器生成的 \(\alpha\) 穩定器 \(G_{\alpha}\)

\[ G_{\alpha} = \langle t_{i}st_{i^{s}}^{-1} | i \in \alpha^{G}, s \in S \rangle \]

其中 \(t_{i}\) 被定義為將 \(\alpha\) 移到 \(i\)\(G\) 中一個元素,即 \(i\) 的陪集(coset)代表。

強生成集

定義 \(G\)(base)\(B\) 為一個序列 \(B = (b_{1},\cdots,b_{k}) \subseteq \Omega\), 因此逐點穩定器(pointwise stabilizer)\(G_{b_{1},\cdots, b_{k}}\) 是可以忽略不計的。

定義 \(G\) 相對於 \(B\)強生成集(SGS,strong generating set)是一個對於每個 \(i\)\(\langle S \cap G^{(i)}\rangle = G(i)\) 的集合 \(S \subseteq G\),其中 \(G^{(i)} := G_{b1,\cdots,bi}\),\(G^{(0)} := G\)

過程

用來計算排列羣(permutation group)\(G\) 階數的 Schreier–Sims 算法一般有三種。它們同樣也可以用來計算 \(G\) 的基和強生成集。

基礎 Schreier–Sims 算法

  1. 如果 \(G\) 為不平凡(non-trivial)羣,選擇一個尚未被選擇的點 \(b\in \Omega\)
  2. 計算根為 \(b\) 的 Schreier 樹,得到 \(|b^{G}|\)
  3. 使用 Schreier 引理找到 \(G_{b}\) 的生成器。
  4. \(G_{b}\) 遞歸使此算法,到找到 \(|G_{b}|\) 為止。

當算法結束時,\((b_{1},\cdots,b_{m})\) 是一個基,找到的所有生成器的並集是一個強生成集。

基礎 Schreier–Sims 算法的運行時間是 指數級 的,但可以被優化成更高效的算法。

增量 Schreier–Sims 算法

增量 Schreier–Sims 算法 是常被用來構建強生成集的快速算法。

如果有一個羣 \(G\) 的強生成集,因為已經得到了所有 \(G^{(i)}\) 穩定器的生成器,那麼很容易得到 \(G\) 的階。部分基(partial base)\(B = [b_{1},\cdots, b_{k}]\) 和部分強生成集 \(S\) 是集合 \(B \in \Omega\) 和集合 \(S \subseteq G\),使得 \(S\) 的任何元素都不能固定 \(B\) 的每個元素。

增量 Schreier–Sims 算法可以將任意部分基和部分強生成集轉化為基和強生成集。

定義 \(T_{i+1}\) 為通過 Schreier 樹 \(G^{(i)}\)\(G^{(i+1)}\) 的陪集的作用。

  1. 如果 \(S = {}\),返回 \(B, S\);
  2. 非空部分基 \(B = (b_{1},\cdots, b_{k}]\)。部分強生成集 \(S\)。集 \(C:= [b_{2},\cdots,b_{k}]\),\(T := S \cap G_{b1}\),並遞歸地應用於輸入 \(C, T\),以將它們修改為 \(H = \langle T \rangle\) 的基和強生成集。
  3. \(B := B \cup C\),\(S := S \cap T\)。用篩選算法在 \(H \leqslant G_{b_{1}}\) 中進行成員資格測試(Membership testing,檢查集合(列表、集合、字典等)是否包含特定元素)。對 \(G_{b_{1}}\) 測試每個 Schreier 生成器 \(s\) 以查看 \(s \in H\)。如果都在 \(H\) 中,那麼有 \(H = G_{b_{1}}\), 返回 \(B,S\)。否則到步驟 4。
  4. 否則有一個 Schreier 生成器 \(s \in G_{b_{1}}\)\(s \notin H\)。設 \(S := S \cup {s}\)。如果 \(s\) 固定了 \(B\) 的所有點,將一個由 \(s\) 移動的 \(\Omega\) 點附加到 \(B\)。回到步驟 2。

當算法結束時,\(B\) 為基,\(S\) 是大小為 \(O(n^{2}\log n)\) 的強生成集。

增量 Schreier–Sims 算法的運行時間為 \(O(n^{8}\log^{3}n)\),即 \(n\) 的多項式。\(t\) 個生成器構建 Schreier 樹需要 \(O(n^{2}+nt)\),或對於 \(t>n\)\(O(nt)\)。因為已經用 \(O(n^{2}\log n)\) 限制了 Schreier 生成器 \(t\) 的數量,所以每個篩選過程都可以在 \(nO(n(n^{2}\log n)) = O(n^{4}\log n)\) 中完成。

實現

以下為基礎 Schreier–Sims 算法的參考代碼。

參考代碼
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#include <iostream>
using namespace std;

const int maxn = 50;     // Maximum size of omega = {1, ,n}
const int maxr = 10000;  // Maximum number of generators

class Permutation {  // interface for permutations
 public:
  int p[maxn];  // the images of the points 0..   maxn-1

  Permutation() { n = maxn; };  // constructors

  Permutation(int m) { n = m; };

  Permutation(int m, char c) {
    n = m;
    switch (c) {
      case 'i':
        for (int i = 0; i < n; i++) p[i] = i;
        break;  // identity
      default:
        for (int i = 0; i < n; i++) p[i] = -1;
        break;  // undefined
    }
  }

  Permutation operator*(Permutation param) const {  // multiplication
    Permutation result(n);
    for (int i = 0; i < n; i++) result.p[i] = param.p[p[i]];
    return (result);
  }

  void operator*=(Permutation param) {  // direct multiplication
    for (int i = 0; i < n; i++) p[i] = param.p[p[i]];
  }

  Permutation inverse() const {  // inverse
    Permutation result(n);
    for (int i = 0; i < n; i++) result.p[p[i]] = i;
    return (result);
  }

  bool isdefined() const { return (p[0] > -1); }  // if it is     defined

  bool isidentity() const {  // if it is the     identity
    for (int i = 0; i < n; i++)
      if (p[i] != i) return false;
    return true;
  }

  bool operator==(Permutation param) const {  // comparison
    for (int i = 0; i < n; i++)
      if (param.p[i] != p[i]) return false;
    return true;
  }

  void input() {
    for (int i = 0; i < n; i++) {
      cin >> p[i];
      p[i]--;
    }
  }  // input

  void output() const {
    for (int i = 0; i < n; i++) cout << p[i] + 1 << " ";
    cout << endl;
  }  // output

  void setn(int m) { n = m; }

 private:
  int n;  // size of omega = {1, ,n}
};

int n;                                   // size of omega = {1, ,n}
int r;                                   // number of generators
Permutation* g = new Permutation[maxr];  // the generators
int nr;
Permutation* newg = new Permutation[maxr];
int cosreps;  // number of    (= size of orbit of alpha)
Permutation* cosrep =
    new Permutation[maxn];  // coset    representatives (to store the output of
                            // SchreierTree)
Permutation undefined(maxn, 'u');

/****** ScheierTree ******/
void ScheierTree(
    int alpha) {  // depth first search to determine the orbit of     alpha
  static Permutation gen(
      n, 'i');    // group element moving original alpha to actual    alpha
  static int ag;  // image of actual alpha under generator g[j]
  cosrep[alpha] = gen;
  cosreps++;
  for (int j = 0; j < r; j++) {
    ag = g[j].p[alpha];
    if (!cosrep[ag].isdefined()) {
      gen *= g[j];
      ScheierTree(ag);
      gen *= g[j].inverse();
    }
  }
}

void SchreierSims() {
  int alpha = 0;
  Permutation sg;
  cout << "THE ORDER OF THE GROUP:\n";
  do {
    cosreps = 0;
    for (int i = 0; i < n; i++) cosrep[i] = undefined;
    // get the coset representatives for G(alpha)
    ScheierTree(alpha);
    // schreier lemma loop to get the schreier generators
    nr = 0;
    for (int i = 0; i < n; i++) {
      if (cosrep[i].isdefined())
        for (int j = 0; j < r; j++) {
          // calculate the (schreier) generators
          sg = cosrep[i] * g[j] * cosrep[g[j].p[i]].inverse();
          bool alreadyhavethis = sg.isidentity();
          for (int k = 0; k < nr; k++)
            if (sg == newg[k]) alreadyhavethis = true;
          if (!alreadyhavethis) newg[nr++] = sg;
        }
    }
    cout << cosreps << flush;
    if (cosreps > 1) cout << "*";
    r = 0;
    for (int j = 0; j < nr; j++) {
      g[r++] = newg[j];
    }
    alpha++;
  } while (cosreps > 1);
  cout << endl;
}

int main() {
  cout << "n ( Size of Omega = {1..n} ) ? ";
  cin >> n;
  for (int j = 0; j < n; j++) {
    g[j].setn(n);
    newg[j].setn(n);
  }
  undefined.setn(n);

  cout << "How many group generators ? ";
  cin >> r;
  for (int j = 0; j < r; j++) g[j].input();

  SchreierSims();
  delete[] g;
  delete[] newg;
  delete[] cosrep;
  return 0;
}

例題

Grand Prix of Yekaterinburg 2015 Problem H Heimdall

海姆達爾——阿斯加德最偉大的兒子之一,眾神和世界之樹的守護者。自古以來古他的主要職責就是守衞阿斯嘉德的入口——一座世界之間的橋樑。現存唯一古老的技術是將一定數量的橋樑結合起來,創造出一座穿越中間世界的橋樑。例如:如果第一座橋將物質從世界 A 傳輸到世界 B,第二座橋——從 B 到 C,那麼它們的組合可以直接將物質從世界 A 傳輸到世界 C. 而且,這個古老的技術甚至可以讓你自己結合一座橋。海姆達爾想知道——使用他所知道的橋樑以及它們的組合,可以創造出多少不同的橋樑。輸入兩個整數 \(R\),\(N\) 分別是海姆達爾發現的橋樑總數和宇宙中的世界數(\(1 \leqslant N \leqslant 15\),\(1 \leqslant R \leqslant 1000\))。接下來的 R 行包含這些橋的信息。每個橋由 \(N\) 個整數 \(a_{1}, a_{2},\cdots a_{n}\) 組成。其中 \(a_{i}\) 表示物質可以通過當前的橋樑轉移到世界 \(i\)。如果當前的橋不影響那些世界,\(a_{i} = i\)。請輸出一個可以通過古老技術建造的不同橋樑的總數。

參考資料與拓展閲讀

  1. Schreier–Sims algorithm - Wikipedia
  2. Knuth, Donald E,Efficient representation of perm groups, Combinatorica 11 (1991), no. 1, 33–43.
  3. Ákos Seress,Permutation Group Algorithms, Cambridge University Press
  4. Sims, Charles C,Computational methods in the study of permutation groups, Computational Problems in Abstract Algebra, pp. 169–183, Pergamon, Oxford, 1970.
  5. Martin Jaggi,Implementations of 3 Types of the Schreier-Sims Algorithm, MAS334 - Mathematics Computing Project, 2005
  6. The Schreier-Sims algorithm for finite permutation groups
  7. Henrik B¨a¨arnhielm,The Schreier-Sims algorithm for matrix groups