むかしのコード（凸多角形の三角形分割とDSS）
今回の記事の目的は、むかし書いた恥ずかしいコードをさらすことで、人見知りを克服しようというものです。参考にしていたサイトがなくなっていたので覚書として公開させて下さい。よろしくお願いします。
むかし、Davenport-Schinzel 列を自習していたときのコードです。「自然数nが与えられたとき、凸n角形の三角形分割をすべて列挙せよ」とても読みにくいんですけど、読めば何となく分かると思います。Mac でコンパイルできるよう少し編集しましたが、基本そのままです。汚いコードが好きな方はどうぞご覧ください。なんと驚きの180行！なげぇ。
関数 make_tree が三角形分割の家系図を作り、traversal で葉っぱの三角形分割を描画しつつDSSをファイルに出力します。下記のコードは、8角形に対するもので、長さが13の(7, 3)-DS 列を出力します*1。
また、16角形で実行しようとしたらセグメンテーションフォールトで落ちたってメモが残してありました。きっと「未来の自分がなんとかしろよ」的な意味合いで残したのでしょう。うざいですねぇ。free すら用意しないくせにですよねぇ。
参考情報についてです。凸多角形の三角形分割の列挙アルゴリズムの文献は不明です。以前は下記コードの最初のコメントのリンクの解説を参考にしました。１時間位でサクッと書き上げられるぐらい良くまとまっていたのに。少し残念。DSSはDavenport?Schinzel Sequences and their Geometric Applications（これから買うならペーパーバックでもよろしいかと）。それから、それらの個数は木の個数と同じなのでグラフ理論入門を見ています。関連構造は数え上げ組合せ論入門 (日評数学選書) にあります。Yes, Catalan 数。
図形と文字列（と木構造など）の関係性がかなり面白い大好きな例題です。構造のコード化って、すっごくロマンチックですよね！
/* * Enumerating triangurations for a convex n-gon, and * their induced (n-1, 3) Davenport-Schinzel sequences. * Ref. http://www.cs.mcgill.ca/~arajwa/compgeom/ */ #include <glut.h> #include <fstream> #include <cmath> using namespace std; enum { X = 0, Y = 1, unconnected = 0, connected = 1, Kgon = 8 }; typedef struct tri_tree { int n; int degree; unsigned char is_leaf; unsigned char adj_mat[Kgon][Kgon]; struct tri_tree **sons; } T_TREE; /*** globals ***/ double angular[Kgon][2]; T_TREE *root = NULL; int ID = 0; int num = 0; void make_tree(T_TREE *parent); void set_angular(void); void traversal(T_TREE *t_ptr, ofstream& fout); void init(void) { int i, j; root = (T_TREE*) malloc(sizeof(T_TREE)); root->n = 3; root->degree = 2; if (Kgon == 3) root->is_leaf = 1; else root->is_leaf = 0; for (i=0; i<3; i++) { for (j=0; j<3; j++) { root->adj_mat[i][j] = unconnected; } } root->adj_mat[0][1] = root->adj_mat[1][0] = connected; root->adj_mat[1][2] = root->adj_mat[2][1] = connected; root->adj_mat[2][0] = root->adj_mat[0][2] = connected; root->sons = (T_TREE**) malloc(sizeof(T_TREE*) * 2); root->sons[0] = root->sons[1] = NULL; set_angular(); make_tree(root); } void make_tree(T_TREE *parent) { int i, j, k, n, deg; T_TREE *t_ptr; if ((n = parent->n) == Kgon) return; for (i=0, deg=0; i<n-1; i++) { if (parent->adj_mat[n-1][i]) { t_ptr = (T_TREE*) malloc(sizeof(T_TREE)); t_ptr->n = n+1; /* init adj matrix */ for (j=0; j<=n; j++) for (k=0; k<=n; k++) t_ptr->adj_mat[j][k] = unconnected; /* copying diagnals not connecting parent of v_{n-1} */ for (j=0; j<n-1; j++) { for (k=0; k<n-1; k++) { if (parent->adj_mat[j][k] == connected) t_ptr->adj_mat[j][k] = connected; } } /* if vertex v_k(0<=k<=i) connecting v_{n-1} exists, * then add (v_k, v_n). */ for (j=0, k=0; j<=i; j++) { if (parent->adj_mat[n-1][j] == connected) { t_ptr->adj_mat[n][j] = t_ptr->adj_mat[j][n] = connected; k++; } } /* if vertex v_k'(i<=k'<=n-2) connecting v_{n-1} exists, * then add (v_k', v_{n-1}). */ for (j=i; j<n; j++) { if (parent->adj_mat[n-1][j] == connected) t_ptr->adj_mat[n-1][j] = t_ptr->adj_mat[j][n-1] = connected; } t_ptr->adj_mat[n-1][n] = t_ptr->adj_mat[n][n-1] = connected; t_ptr->degree = k+1; if (n+1 == Kgon) { t_ptr->is_leaf = 1; } else { t_ptr->is_leaf = 0; t_ptr->sons = (T_TREE**) malloc(sizeof(T_TREE*) * (k+1)); for (j=0; j<=k; j++) { t_ptr->sons[j] = NULL; } } parent->sons[deg++] = t_ptr; } } if (n+1 != Kgon) { for (i=0; i<deg; i++) make_tree(parent->sons[i]); } } void set_angular(void) { int i; const double d = 2 * M_PI / Kgon; double ang; for (i = 0, ang = 0.; i < Kgon; i++, ang += d) { angular[i][X] = cos(ang); angular[i][Y] = sin(ang); } } void traversal(T_TREE *t_ptr, ofstream& fout) { int i, j; if (t_ptr == NULL) return; if (t_ptr->is_leaf == 1) { num++; for (i=0; i<Kgon-1; i++) { for (j=i+1; j<Kgon; j++) { if (t_ptr->adj_mat[i][j] == connected) { glColor3f(0.0, 0.0, 0.0); glBegin(GL_LINES); glVertex2d(angular[i][X]+2.5*(ID%12),angular[i][Y]+2.5*(ID/12)); glVertex2d(angular[j][X]+2.5*(ID%12),angular[j][Y]+2.5*(ID/12)); glEnd(); } } } ID++; for (i=1; i<Kgon; i++) { for (j=i-1; j>-1; j--) { if (t_ptr->adj_mat[i][j] == connected) { fout.width(2); fout << j; } } } fout << endl; return; } for (i=0; i<t_ptr->degree; i++) { traversal(t_ptr->sons[i], fout); } } void display(void) { glClear(GL_COLOR_BUFFER_BIT); ofstream fout("dss.out"); ID=0; traversal(root, fout); fout.close(); glutSwapBuffers(); } void resize(int w, int h) { glViewport(0, 0, w, h); glLoadIdentity(); gluOrtho2D(-1.5, 30, -1.5, 30); } int main(int argc, char* argv[]) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE); glutInitWindowSize(600, 600); glutInitWindowPosition(250, 0); glutCreateWindow("Enumerate All Triangulations of a convex n-gon."); glClearColor(1.0, 1.0, 1.0, 1.0); glEnable(GL_LINE_SMOOTH); glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glHint(GL_LINE_SMOOTH_HINT, GL_NICEST); glutReshapeFunc(resize); glutDisplayFunc(display); init(); glutMainLoop(); return 0; }
- 実行結果
下の画像は、家系図を作った後、深さ優先順で出力している例です。左下から右上に向かってラスター順に並んでいます。
[
:W555]
以下にDSSの例を示しています。まず、多角形の最右の頂点番号が0でそこから反時計回りの順で番号が振られています。「DSS→三角形分割」の見方は、DSS を単調減少部分列で分割してみます。それを順番に、1, 2, ..., n-1と番号をつけます。この番号と各頂点番号が対応づいていて、部分列はどの頂点と接続しているかを表しています。例えば、上記の画像の左下の三角形分割とDSS s = 0 10 20 30 40 50 60 は相関関係があります。
0102030405060 0102030405650 0102030454060 ... 0123454363210 0123454643210 0123456543210
CC = g++ GLUTIN = /System/Library/Frameworks/GLUT.framework/Headers CFLAGS = -O2 -Wall -I$(GLUTIN) LD = g++ LFLAGS = LIBS = -framework OpenGL -framework GLUT -framework Foundation .SUFFIXES: .SUFFIXES: .o .c PROGS = enum_tri OBJS = $(PROGS:=.o) all: $(PROGS) .c.o: $(CC) $(CFLAGS) -c $< .o: $(LD) $(LFLAGS) $*.o -o $* $(LIBS) .c: $(CC) $(CFLAGS) -c $< $(LD) $(LFLAGS) $(LIBS) $*.o -o $* clean: rm $(PROGS) $(OBJS)
- おしまい
もっと読みやすいコードにしろとか、説明メモ残せとか、家系図見せろとか、フリップグラフ出せとか、いろいろ思います。特に、フリップグラフは描いてみたいですよね。DSS のリストがあれば起こせそうです。描くかどうかは置いておいても、直径とか幅とかいろいろ見てみたい気はしますよね。どっかの文献に記載されているのかもしれませんが。ランダムフリップのカバータイムとかも面白そうですよね。どっかの文献に記載されてそうですが。
*1:楽しみのためあまり丁寧な説明はしたくありません。でちょ補足。列の長さは多角形の辺数nと三角形分割の対角線数n-3の和です。またDSSのアルファベットは全順序関係があり、各アルファベットの最初の出現位置の順序はアルファベットの順序を違えてはいけない。設問は、この2つの条件を満たすDSSをすべて列挙することと等価です。また、以降あまりきちんとした日本語を使っていません。文字列と図形をいろいろ見比べてみて相関関係が分かると、きっと楽しいと思います
