プロファイル

今日はプロファイルについて書きたいと思います。プロファイルとは、プログラムの中の個々の関数がどれくらい呼び出されているか、そしてどれくらい時間を消費しているかの統計を調べることを言います。具体的には、Cでコンパイルするときに、-pgオプションをつけてコンパイルをします。その後gprofプログラムにより、プロファイルを表示することができます。

$ gcc -pg demo.c subgroup.c subset.c list.c cayley.c permutation.c
$ gprof ./a.out gmon.out

例えば次のサンプルは、対称群（置換群）の部分群を求めるプログラムをプロファイルしたものです。このプロファイルからは、is_set0とq_sortにおいて多くの時間が割かれていることがわかります。従って、プログラムの速度を上げるには、この2つの関数を改良することが効果的です。実際にq_sorをCの標準ライブラリのqsortに代え、is_set0（これはセットであるかどうかのチェックをする関数）をノーチェック、つまり常にセットが正常だという返り値を戻すようにしたところ、60秒かかっていたものが14秒に短縮されました。

Each sample counts as 0.01 seconds.
  %   cumulative   self               self     total           
 time   seconds   seconds     calls   s/call   s/call  name    
 44.64     25.84    25.84 381729969     0.00     0.00  is_set0
 40.11     49.05    23.21     29581     0.00     0.00  q_sort
  4.27     51.52     2.47 190741601     0.00     0.00  set_cmp
  2.66     53.05     1.54     59320     0.00     0.00  l_lookup_with_cmp
  2.28     54.38     1.32 381761640     0.00     0.00  fact
  2.27     55.69     1.31 381670649     0.00     0.00  is_set
  1.02     56.28     0.59     23761     0.00     0.00  l_lookup
  0.92     56.81     0.53     14032     0.00     0.00  l_append
  0.76     57.25     0.44     29581     0.00     0.00  uniq
  0.59     57.59     0.34     29581     0.00     0.00  set_mul_ss
          …  …  …  …  …  …  …  …  …  …
          …  …  …  …  …  …  …  …  …  …

このようにプロファイルを調べることは実行速度を上げるためには便利です。また、is_set0のようなチェック機能は、堅いプログラム（バグの少ないプログラム、ロブストともいう）では必要ですが、一旦動くことが確認できれば、チェックをはずすことも速度向上につながります。

Symmetric_group レポジトリを更新

対称群の部分群を求めるプログラム（symmetric group）をCで書き直しました。Ruby ではなかなか速度が上がらないので、Cに書き直せば速度が上がると期待したのです。ところが、実際にはさして速度が上がりませんでした。5次対称群の部分群を求めるプログラムでは、最初の段階でRubyが80秒、Cが60秒ぐらいでした。よく調べてみると無駄な箇所があったので、RubyもCも書き換えたところ、Rubyが1分を切り（56秒）、Cは40秒で答えが出るようになりました。

しかしまだ満足はしていません。インターネットで対称群の部分群を求めるプログラムを探してみましたが、アルゴリズムを解説しているようなものは見つかりませんでした。おそらく、さほど新しい問題ではないので、すでに誰かが良いアルゴリズムを見つけているのではないかと思います。特に有限群の研究者にとっては、コンピューターがかなり威力を発揮すると思うので、すでにこの問題の有力なアルゴリズムが既知のものになっているのではないでしょうか。

今日、Githubのリポジトリを更新しました。興味のある人はここを見て下さい。また。良いアルゴリズムをご存知の方がいたら教えていただければ助かります。最後に4次対称群と5次対称群の部分群の個数について記しておきます。

• 4次対称群 => 30個
• 5次対称群 => 156個

クイックソートを書いていてハマった話

今日はクイックソートでハマったことを書きたいと思います。クイックソートはデータを小さい順（または大きい順）に並べるアルゴリズムです。ここでは与えられたデータを「集合」と呼びますが、数学で扱う集合とは異なり、同じ要素を複数含んでいても良いことにします。

クイックソートでは、データの集合Sを3つの部分に分割します。そのためにある基準となる値（ピボットと呼ばれる）xを決めます。xは集合Sの中から取ってくることにします。第一の集合は、xよりも小さい要素の集合S1です。第二は、データxと同じ値を持つ要素の集合S2。第三はxよりも大きい値を持つ要素の集合S3です。S2は同じ値の要素だけなので、すでに整列されています。次の段階は、S1とS3について同じ方法でソートすることです。これらの集合は要素の数がSより小さいので、同じアルゴリズムの有限回の繰り返しで整列が完了します。

このアルゴリズムは完全なのですが、3分割を2分割で済ませられれば、実装が簡単になります。例えば、S1は同じで、S2とS3を一つの集合S4にするのです。 つまりSがxより小さい要素の集合、S4がx以上の要素の集合になります。 ところがこれだとうまくいかない場合があるのです。xの取り方によってはSの全ての要素がx以上であることがあり得ます。その場合S1は空集合、S4はSに等しくなります。ということは、S4の要素の数が減りませんから、アルゴリズムを繰り返すと無限ループに陥ってしまいます。このバグに気づくまで、何時間もかかってしまいました(TT)。

最初の方法では、xがSの要素で、S2が空ではないことが効いています。第二の方法では何らかの方法で、無限ループを回避することが必要です。例えば、S1が空集合だった場合に、S4をS2とS3に分けることにより、第一の方法に帰着させることが考えられます。この方法で書いたCプログラムを以下に載せておきます。何か間違いがあったら、コメントで指摘していただくと助かります。

クイックソート：　配列a[]のインデックスlからhまでの要素を整列する。

void
q_sort(int a[], int l, int h) {
  int swap, l1, h1, s;

  if (l >= h)
    return;
  s = a[(l+h)/2];
  l1 = l;
  h1 = h;
  while (l1 <= h1) {
    for (;a[l1] < s; ++l1)
      ;
    for (;a[h1] >= s;--h1)
      ;
    if (l1 < h1) {
      swap = a[l1];
      a[l1++] = a[h1];
      a[h1--] = swap;
    }
  }
  if (l1 == l) {
    h1 = h;
    while (l1 < h1) {
      for (;a[l1] == s; ++l1)
        ;
      for (;a[h1] > s;--h1)
        ;
      if (l1 < h1) {
        swap = a[l1];
        a[l1++] = a[h1];
        a[h1--] = swap;
      }
    }
  }
  if (h1 == h)
    return;
  q_sort (a, l, h1);
  q_sort (a, l1, h);
}

追記：

投稿後、「プログラム言語 C」のクイックソートプログラムを見たら、発想は同じだが、もっとスマートなプログラムでした。当然ですね。当時の最先端の研究者のプログラムですから。

この方法は、xがS4の最小値であることを利用するものです。はじめに、配列中のxを値に持つ要素を配列の端（プログラム言語Cでは左端）と交換して、残りの部分で集合S1と集合S4に分けます。そのあと、はじめに端に置いた要素と、S1の右端の要素（最初に左端においた場合）を交換、またはS2の左端の要素（最初に右端においた場合）と交換します。あらためて、xより小さい要素の集合をS1、最後に端から中間に移動した要素を除いたx以上の要素の集合をS2とします。

(1)   S1={s|s<x}
(2)   端から中間に移動した要素
(3)   S4={s|s>=x,ただし、端から中間に移動した要素を除く}

これで、S4がSよりも少ない要素数になり、無限ループを避けることができます。

やはり、基本は3つの部分に分けることだったんですね。「プログラム言語C」のクイックソートをそのまま載せるわけにはいかないので、新たにプログラムを書きました。このプログラムでは右端にピボット要素を移動しています。

追記の追記：

申し訳ありません。プログラムにバグがありましたので、修正しました。

• 14行目のwhile文の不等号<を<=に修正しました。
• 15行目と17行目のfor文の条件に、15行目に「&& l1<=h-1」17行目に「&& h1 >= l」を入れました。

static void
q_sort(int a[], int l, int h) {
  int swap, l1, h1, s, mid;

  if (l >= h)
    return;
  mid = (l + h) / 2;
  s = a[mid];
  swap = a[h];
  a[h] = a[mid];
  a[mid] = swap;
  l1 = l;
  h1 = h-1;
  while (l1 <= h1) {
    for (;a[l1] < s && l1<=h-1; ++l1)
      ;
    for (;a[h1] >= s && h1 >= l;--h1)
      ;
    if (l1 < h1) {
      swap = a[l1];
      a[l1++] = a[h1];
      a[h1--] = swap;
    }
  }
  swap = a[h];
  a[h] = a[l1];
  a[l1++] = swap;
  q_sort (a, l, h1);
  q_sort (a, l1, h);
}

Githubにsymmetric_groupをアップロード

2/26追記

本日、プログラムを書き直してレポジトリにアップロードしました。demo.rbも書き換えられていて、下にあるプログラムは新バージョンでは動きません。レポジトリに含まれているdemo.rbを使ってください。変更の主な理由は置換の積の速度をあげることです。元バージョンでは置換の積をその定義どおりに計算しており、計算量が多くなっています。それを新バージョンでは、置換を順に並べ、そのインデックスで置換を表し、インデックスの2次元配列の乗算表を用いることにより、格段に積の計算速度を上げています。これで、5次対称群のすべての部分群を求めるのに、1時間10分＝＞１分20秒に短縮されました。最初のアルゴリズムが悪すぎ。（追記以上）。

GithubにSymmetric_groupというプログラムをアップロードしました。このプログラムはRubyで書かれています。n次対称群（これは置換群と言った方が分かりやすいかもしれません）を計算をしたり、生成部分群を求めたりするプログラムです。

置換群の積の計算は、いくつかの数をたどりながら行うので面倒です。それをコンピューターにやらせれば楽になるのではないかと思い、プログラムを作りました。また、このプログラムは単に掛け算をサポートするだけではなく、群の生成をサポートします。例えば4次の対称群において、2つの要素から部分群を生成したいとします。この部分群は2つの要素を含むような最小の部分群です。与えられた2つの置換から始め、その積や逆元を繰り返し計算します。そしてすでに出てきた要素は捨て、新しい要素は追加し、これ以上新しい要素が出なくなったら、それらの要素の全体が部分群になります。手作業でやるのは、少し大きな部分群では面倒であり、場合によっては事実上不可能です。それをコンピューターで計算すれば楽だし、実用的だと思いました。

プログラムのPGroupクラスのインスタンス生成（initializeメソッド）でこの計算をしていますが、非常にシンプルです。このとき有限集合Aが、AA = Aならば、Aは部分群になるという定理を使っています。

例えば二面体群D4は、（1234）の巡回群と（13）の互換をもとに生成します。これをプログラムで計算すると瞬時に行ってくれます。人間の計算でもそんなに大変ではなのですが、プログラムしたほうがずっと早く答えを見つけられます。また、（1234）の巡回群と（12）の互換で群を生成すると、4次対称群全体になります。これもコンピューターで計算するとあっという間に分かります。demo.rbがそのプログラムで、実行すると次のようなメッセージが出ます。

D4 is a dihedral group. D4 includes rotations and refletions of a square.
It is expressed with a permutation of vertices and it is generated from
two cycles [1,2,3,4] and [1,3],
Let s be [1,2,3,4] and t be [1,3].
They are created with Permutation.new([2,3,4,1]) and
Permutation.new([3,2,1,4]) respectively.

st = PSet.new (s,t)
D4 = PGroup.new (st)

D4 = [[1,2,3,4],[1,3],[1,3][2,4],[1,4][2,3],[1,2][3,4],id,[1,4,3,2],[2,4]]


If you generate a group from [1,2,3,4] and [1,2] instead of the cycles above,
the symmetric group S4 will be generated.
s = Permutation.new([2,3,4,1])
t = Permutation.new([2,1,3,4])
st = PSet.new (s,t)
S4 = PGroup.new (st)

S4 = [[1,2,3,4],[1,2],[1,3][2,4],[1,3,4],[2,3,4],id,[1,4,3,2],[1,4,2,3],
[1,2,4,3],[1,3,2,4],[1,3,4,2],[2,4,3],[1,3,2],[1,4,2],[1,4,3],[1,4][2,3],
[1,4],[3,4],[2,3],[2,4],[1,3],[1,2,3],[1,2,4],[1,2][3,4]]

置換群は、群論の本では例題として良く取り上げられ、面倒な計算が練習問題に良く出されていました。学習者は、その計算を手作業で時間をかけて行っていたのですが、このプログラムを使うと時間短縮になり、よりスムーズに群の学習が進むと思います。

3次元を2次元に変換して表示するプログラムを作成

t2tというアプリケーションをGithubにアップロードしました。

t2tは「Three dimension To Two dimension」を短くしたものです。すなわち、3次元の座標系から2次元の座標系への変換をするアプリケーションです。

例えば地上にある立方体は、空間座標でその各頂点を表すことができます。その物体をドローンから撮影すると物体は2次元のスクリーン上に写されます。元々の3次元の座標をスクリーン上の2次元の座標に変換しているわけです。この計算をするのがt2tです。次のスクリーンショットはt2tの画面です。地上にある立方体をドローンから撮影した画像を示しています。

t2tでは座標変換を使っています。その原理は直交変換を含むアフィン変換と、3次元を2次元に写す相似写像です。詳しいことは省略しますが、「座標変換」などのワードでググるとネット上の情報が得られるでしょう。

t2tではプログラムする必要があります。プログラムは非常に簡単なコマンドをいくつか並べたものです。その詳細はドキュメントで解説されています。 なお、githubのt2tのトップページでは英語の説明が現れますが、日本語の説明が、readme.jp.mdとdoc.jp.mdにあります。上の図の左側のウィンドウに書かれている暗号のような文字列がプログラムです。そして、そのプログラムを実行した結果が右側のウィンドウに表示された図形です。

t2tを使ってみると描かれた図形が歪んで見えることがあります。それは視点（ドローン）と物体の各点を結んでいる直線が平行ではないためです。その直線を平行に近づけるためには視点を遠ざけることが必要です。そこで視点をなるべく遠くにするのですが、今度は物体の像が小さくなってしまいます。そこで適当にスクリーンの画像を大きくすることが必要になってきます。

もう一つの物体を本物らしく見せるための方法は、なるべく多くの物体を配置することです。例えば次の図形は道路に直方体の建物がいくつか並んでいる様子を簡易的に表したものです。 一つ一つの直方体は歪んで見えますが全体としてはそれらしく見えるようになります。

「ブラサガリ算」をコンピュータで解く

　芦ケ原伸之という有名なパズル作家をご存知だろうか。数々のパズル本を出していて、その一つに「超々難問数理パズル」という本がある。この本が出版されたのは2002年で、私が買ったのもおそらくその頃だと思う。購入して、いくつか解いてみたが、当時は仕事も忙しく、やがて書棚に置きっぱなしになってしまった。この本の裏表紙に「この一冊でたっぷり10年は楽しめます！」と書かれているが、私の場合は10年以上（本を開かなかったので）楽しまなかったことになる。退職して時間に余裕ができたので、あらためて取り組んでみた。タイトルの通り、かなり難しい。今回はその中から9番目の問題「ブラサガリ算」という問題について書きたいと思う。

　この問題は逆三角形に並べた数字で、隣り合う2数の差が下に来るという単純なルールが成り立つものを見つけるのである。例えば3個の数字（1,2,3の3つの数）で逆三角形をつくる場合は簡単で、下のような図になる。なお、答えは一種類ではない。

      2 3
       1

　次に6個の数字にすると少し難しくなるが、しばらく試行錯誤を繰り返せばできるだろう。これも答えは一通りではない。

     4 6 1
      2 5
       3

　さて、問題はこれを5段組15個の数字（1から15の整数）で三角形を作れ、というものだった。

　 まず、手作業でやってみたがなかなかできない。大きな数、例えば15や14が上の方に来なければいけないということはわかったが、それでも場合が多くて答えが見つからなかった。それならコンピューターに計算をさせて全数テストをしてしまえ、ということになった。ところが、このプログラムを考えるのも結構難しかったのである。できあがったプログラムはここに記しておくので、興味のある方は、プログラムを走らせて答を見てほしい。

※　以下の枠で囲まれた文は、最初に投稿したときのものである。その後、rubyのArray（配列）オブジェクトに標準でpermutationメソッドがあるのに気がついて、プログラムを書き直した。

　プログラムの中に、n個からr個をとってできる順列をひとつずつ取り出すというメソッドがある。具体的にはプログラムを見てもらいたいが、再起手続きとyield命令を使って、ひとつひとつの順列をブロックで処理できるようなメソッドになっている。この計算は一般的なものなので、すでにライブラリがあるのではないかと思う。もしこれをお読みの方で、そのようなライブラリーを知ってる方がいらっしゃったら、コメントで教えていただければ幸いです。

　プログラムを走らせると1秒ほどで答えが出る。答えは対象を除き、一通りである。ついでに6段の場合も試してみた。2分4秒かかって、条件を満たす並べ方が無いことが分かった。実は、パズル本の解答編にも6段には答えがないと書いてあった。

@n = 5
@max = @n*(@n+1)/2
@a = []
1.upto(@max) do |i|
  @a << i
end
@b = []
0.upto(@n-1) do |i|
  @b[i] = []
end
@c = []

def success?
  0.upto(@n-1) do |i|
    if @c[i] != true
      return false
    end
  end
  return true
end

def check x
  if @c[x] == true
    return false
  end
  @c[x] = true
end

def try
  @c = []
  0.upto(@n-1) do |i|
    unless @b[@n-1][i].instance_of? (Integer)
      raise "@b[#{@n-1}][#{i}] is not an integer.\n"
    end
    if @b[@n-1][i] < 1 || @b[@n-1][i] > @max
      raise "@b[#{@n-1}][#{i}]=#{@b[@n-1][i]} is too big or too small.\n"
    end
    if ! check(@b[@n-1][i])
      return false
    end
  end

  (@n-2).downto(0) do |i|
    0.upto(i) do |j|
      @b[i][j]=(@b[i+1][j+1]-@b[i+1][j]).abs
      if ! check(@b[i][j])
        return false
      end
    end
  end
  true
end

@a.permutation(@n) do |p|
  @b[@n-1] = p
  if try()
    (@n-1).downto (0) do |i|
      p @b[i]
    end
  end
end

音声入力、はじめました

　最近音声入力にハマっている。既に知っている人からは今更何を言ってるんだと思われるかもしれない。しかし自分は最近になって音声入力を始めたので、仕方がないのだ。

　音声入力は二つあって、一つはコンピューターの音声入力で、Chromeを使って入力する（Firefoxなどの他のブラウザでは音声入力はできない）。だから初めにGoogle Chromeをインストールしておく必要がある。Chromeを開いたらドキュメント（グーグル・ドキュメント）を開く。「ツール」メニューから音声入力を選択する。言語は日本語がデフォルトだが、英語に設定すれば英語の音声入力もできるし、他の言語も可能である。昔は口述筆記と言うのがあったが、それには秘書を雇うようなお金が必要だった。今は秘書の代わりをコンピューターが務めてくれるわけだ。

　もうひとつの音声入力はスマホを使うものだ。私が持っているのはAndroidスマホである。入力キーボードはデフォルトで日本語になっているが、それを「Google音声入力」に変えると音声入力が可能になる。これも変換効率は非常に良くて、ストレスなく入力をすることができる。それにキーボードから入力するよりも明らかに早く入力することができる。ただ誤変換が多少あるので、その場合にはキーボードを切り替えて訂正しなければならない。

　最初は慣れなくても、しばらく続けていると慣れてくる。そして、入力が速くなってくる。音声入力に慣れてくると、キーボードからの入力に戻ることはもうできない。楽ちんすぎるのだ。コンピューターが進歩すると、人間が堕落するような気がする。これからはキーボード入力を全然しない人が出てくるかもしれない。