M.Hiroi's Home Page
http://www.geocities.jp/m_hiroi/

Memorandum

プログラミングに関する覚え書や四方山話です。
[ Home | 2018年 1月, 5月, 6月 ]

2018 年 6 月

6月24日

●有限体

数学の世界では、加減乗除ができる数の集合のことを「体」といいます。特に、数が n 個しかないものを「有限体 (finite field)」または「ガロア体 (Galois field)」といいます。ここでは GF(n) と表記することにします。n が素数のとき、その剰余の集合 {0, 1, ..., n - 1} は有限体になります。簡単な例として GF(2) と GF(3) の加法と乗法を示します。

  GF(2)

  + |  0  1    * |  0  1
 ---+-------  ---+-------
  0 |  0  1    0 |  0  0
  1 |  1  0    1 |  0  1

  GF(3)

  + |  0  1  2    * |  0  1  2
 ---+----------  ---+----------
  0 |  0  1  2    0 |  0  0  0
  1 |  1  2  0    1 |  0  1  2
  2 |  2  0  1    2 |  0  2  1

mod n の世界では、加法と乗法は簡単に計算できます。加算または乗算したあと mod n を計算するだけです。減算は負数 (加法の逆元) を、徐算は逆数 (乗法の逆元) を考えます。x + y = 0 を満たす y が x の負数になるので、y = n - x であることがわかります。0 以外の数の逆数が存在すれば徐算が可能になります。乗法の表を見れば、GF(2) は 1 * 1 = 1, GF(3) は 1 * 1 = 1, 2 * 2 = 1 となるので、徐算が成り立つことがわかります。

n が素数ではない場合、逆数が存在しない数があります。たとえば、mod 4 の加法と乗法は次のようになります。

  mod 4 の加法と乗法

  + |  0  1  2  3     * |  0  1  2  3   
  --+-------------   ---+-------------  
  0 |  0  1  2  3     0 |  0  0  0  0
  1 |  1  2  3  0     1 |  0  1  2  3
  2 |  2  3  0  1     2 |  0  2  0  2
  3 |  3  0  1  2     3 |  0  3  2  1

mod 4 の場合、加算、乗算、減算は成り立ちますが、徐算で 2 の逆数が存在しないことが乗法の表からわかります。つまり、有限体にはならないわけです。

●拡大体

4 の剰余 (0, 1, 2, 3) は GF(4) になりませんが、GF(2) に新しい数を付け加えることで、数が 4 つある有限体を作ることができます。新しい数を付加して体を拡張することを「体の拡大」といい、拡大された体を「拡大体」といいます。たとえば、実数の世界では代数方程式 x2 + 1 = 0 に解はありませんが、虚数 i を導入した複素数の世界では解くことができますね。複素数の世界は実数の拡大体になるわけです。

GF(2) を拡大する場合、GF(2) 上の代数方程式を考えます。2 次方程式を列挙すると次のようになります。

  1. x2 = 0
  2. x2 + 1 = 0
  3. x2 + x = 0
  4. x2 + x + 1 = 0

この中で解を持たないのは式 4 です。『x の関数 f(x) において、f(a) = 0 ならば (x - a) で割り切れる』という性質を「因数定理」といいます。式 1, 2, 3 を因数分解すると次のようになります。

  1. x2 = (x + 0)(x + 0)
  2. x2 + 1 = (x + 1)(x + 1)
  3. x2 + x = (x + 0)(x + 1)

式 4 は GF(2) で因数分解することはできません。これを「既約多項式」といいます。既約多項式は多項式の素数みたいなもので、1 次多項式 x, x + 1 で割り切れることはありません。この式の解の一つを a としましょう。a2 + a + 1 = 0 が成り立つので、a2 = a + 1 であることがわかります。次に、a の累乗を計算します。

a1 = a
a2 = a + 1
a3 = a * a2 = a * (a + 1) = a + 1 + a = 1
a4 = a1+3 = a
a5 = a2+3 = a + 1
a6 = a3+3 = 1

a, a + 1, 1 が順番に現れています。これを多項式で表すと 1, x, x + 1 になり、1 次多項式がすべて現れていることがわかります。このように、累乗を計算すると n 次未満の多項式がすべて現れる既約多項式のことを「原始多項式」といいます。これに 0 を加えた 4 つの要素 {0, 1, x, x + 1} を数として加法と乗法を考えてみましょう。基本的には式を計算して、それを既約多項式 x2 + x + 1 で割った剰余が値になります。下図を見てください。

 +  |  0   1   x  x+1     *  |  0   1   x  x+1
----+-----------------   ----+------------------
 0  |  0   1   x  x+1     0  |  0   0   0   0
 1  |  1   0  x+1  x      1  |  0   1   x  x+1
 x  |  x  x+1  0   1      x  |  0   x  x+1  1
x+1 | x+1  x   1   0     x+1 |  0  x+1  1   x  

加法の場合は特に簡単で、係数ごとに加算して mod 2 をとる、つまり係数ごとの排他的論理和になります。乗算の場合は累乗の形式に変換して計算すると簡単です。つまり、am * an = an+m = a(n+m) mod 3 になります。a2 = a + 1 の関係を使って式を簡約してもかまいません。

減法の場合、加法表より負数が自分自身であることがわかるので、加法と同様に係数ごとの排他的論理和になります。徐算の場合、乗法表より 1, x, x + 1 の逆数がすべて存在するので、徐算が成り立つことがわかります。つまり、{0, 1, x, x + 1} は有限体になるわけです。多項式の係数をビットで表すと、加法と乗法は次のようになります。

 +  |  00  01  10  11     *  |  00  01  10  11
----+-----------------   ----+------------------
 00 |  00  01  10  11     00 |  00  00  00  00
 01 |  01  00  11  10     01 |  00  01  10  11
 10 |  10  11  00  01     10 |  00  10  11  01
 11 |  11  10  01  00     11 |  00  11  01  10  

p を素数とし、n を自然数とすると、要素数が pn の有限体は必ず存在することが知られています。これを GF(pn) と表記します。GF(2n) はコンピュータとの相性が良いので、エラー検出、乱数の生成、暗号など多くの分野で応用されています。

ご参考までに、GF(23) と GF(32) の加法表と乗法表を示します。

GF(23)
原始多項式 f(x) = x3 + x + 1 
f(a) = 0 となる要素 a を考える
a1 = a (010)
a2 = a2 (100)
a3 = a + 1 (011)
a4 = a2 + a (110)
a5 = a2 + a + 1 (111)
a6 = a2 + 1 (101)
a7 = 1 (001)
GF(23) の加法表
  
 +  | 000 001 010 011 100 101 110 111 
----+---------------------------------
000 | 000 001 010 011 100 101 110 111 
001 | 001 000 011 010 101 100 111 110 
010 | 010 011 000 001 110 111 100 101 
011 | 011 010 001 000 111 110 101 100 
100 | 100 101 110 111 000 001 010 011 
101 | 101 100 111 110 001 000 011 010 
110 | 110 111 100 101 010 011 000 001 
111 | 111 110 101 100 011 010 001 000 


GF(23) の乗法表, (n) の n は指数

              (7) (1) (3) (2) (6) (4) (5)
   *    | 000 001 010 011 100 101 110 111
--------+----------------------------------
000     | 000 000 000 000 000 000 000 000 
001 (7) | 000 001 010 011 100 101 110 111 
010 (1) | 000 010 100 110 011 001 111 101 
011 (3) | 000 011 110 101 111 100 001 010 
100 (2) | 000 100 011 111 110 010 101 001 
101 (6) | 000 101 001 100 010 111 011 110 
110 (4) | 000 110 111 001 101 011 010 100 
111 (5) | 000 111 101 010 001 110 100 011 
GF(32)
原始多項式 f(x) = x2 + x + 2
f(a) = 0 となる要素 a を考える
a1 = a                                (1 0)
a2 = -a - 2 = 2a + 1                  (2 1)
a3 = a(2a + 1) = 4a + 2 + a = 2a + 2  (2 2)
a4 = a(2a + 2) = 4a + 2 + 2a = 2      (0 2)
a5 = 2a                               (2 0)
a6 = 2a * a = 2(2a + 1) = a + 2       (1 2)
a7 = a(a + 2) = 2a + 1 + 2a = a + 1   (1 1)
a8 = a(a + 1) = 2a + 1 + a = 1        (0 1)
GF(32) の加法表

  +   | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2) 
------+-------------------------------------------------------
(0 0) | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2) 
(0 1) | (0 1) (0 2) (0 0) (1 1) (1 2) (1 0) (2 1) (2 2) (2 0) 
(0 2) | (0 2) (0 0) (0 1) (1 2) (1 0) (1 1) (2 2) (2 0) (2 1) 
(1 0) | (1 0) (1 1) (1 2) (2 0) (2 1) (2 2) (0 0) (0 1) (0 2) 
(1 1) | (1 1) (1 2) (1 0) (2 1) (2 2) (2 0) (0 1) (0 2) (0 0) 
(1 2) | (1 2) (1 0) (1 1) (2 2) (2 0) (2 1) (0 2) (0 0) (0 1) 
(2 0) | (2 0) (2 1) (2 2) (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) 
(2 1) | (2 1) (2 2) (2 0) (0 1) (0 2) (0 0) (1 1) (1 2) (1 0) 
(2 2) | (2 2) (2 0) (2 1) (0 2) (0 0) (0 1) (1 2) (1 0) (1 1)   


GF(32) の乗法表, (n) の n は指数

                   (8)   (4)   (1)   (7)   (6)   (5)   (2)   (3)
  *       | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2)
----------+-------------------------------------------------------
(0 0)     | (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) 
(0 1) (8) | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2) 
(0 2) (4) | (0 0) (0 2) (0 1) (2 0) (2 2) (2 1) (1 0) (1 2) (1 1) 
(1 0) (1) | (0 0) (1 0) (2 0) (2 1) (0 1) (1 1) (1 2) (2 2) (0 2) 
(1 1) (7) | (0 0) (1 1) (2 2) (0 1) (1 2) (2 0) (0 2) (1 0) (2 1) 
(1 2) (6) | (0 0) (1 2) (2 1) (1 1) (2 0) (0 2) (2 2) (0 1) (1 0) 
(2 0) (5) | (0 0) (2 0) (1 0) (1 2) (0 2) (2 2) (2 1) (1 1) (0 1) 
(2 1) (2) | (0 0) (2 1) (1 2) (2 2) (1 0) (0 1) (1 1) (0 2) (2 0) 
(2 2) (3) | (0 0) (2 2) (1 1) (0 2) (2 1) (1 0) (0 1) (2 0) (1 2) 

●参考文献・URL


6月16日

●不定方程式と拡張ユークリッドの互除法

変数の個数が式の本数よりも多い方程式を「不定方程式」といいます。一次不定方程式 ax + by = c (c != 0) は、a, b, c が整数で、c が a と b の最大公約数の倍数であるとき、整数解を持つことが証明されています。これを「ベズーの等式」といいます。ベズーの等式は「拡張ユークリッドの互除法」を使って簡単に解くことができます。

拡張といっても原理はユークリッドの互除法とまったく同じです。gcd(a, b) を求めるときに一組の解 (x, y) もいっしょに見つける、というものです。ポイントは数を x = xa * a + xb * b の形で表すところです。次の図を見てください。

  x = xa * a + xb * b, y = ya * a + yb * b とおき、
  x, y に対してユークリッドの互除法を適用する
  (初期値は xa = 1, xb = 0, ya = 0, yb = 1)

  x と y の商を q, 余りを z とすると

  z = x - q * y = xa * a + xb * b - q * (ya * a + yb * b)
    = (xa - q * ya) * a + (xb - q * yb) * b
    = za * a + zb * b

  となる
  
  次は y と z で同様の処理を繰り返す
  具体的には変数の値を

  x, xa, xb = y, ya, yb
  y, ya, yb = z, za, za

  のように更新して処理を繰り返す
  
  y が 0 になったとき、x が a と b の最大公約数で (xa, xb) が一つの解となる

具体的な例を示しましょう。4321 * x + 1234 * y = gcd(4321, 1234) を解いてみます。

  (x, xa, xb)     (y, ya, yb)      q   (z, za, zb)
  ------------------------------------------------------
  (4321, 1, 0)    (1234, 0, 1)     3   (619, 1, -3)
  (1234, 0, 1)    (619, 1, -3)     1   (615, -1, 4)
  (619, 1, -3)    (615, -1, 4)     1   (4, 2, -7)
  (615, -1, 4)    (4, 2, -7)       153 (3, -307, 1075)
  (4, 2, -7)      (3, -307, 1075)  1   (1, 309, -1082)
  (3, -307, 1075) (1, 309, -1082)  3   (0, -1234, 4321)
  (1, 309, -1082) (0, -1234, 4321)

  4321 * 309 + 1234 * (-1082) = 1 になる

一般解は次のように求めることができます。

   4321 * x   + 1234 * y       = 1
-  4321 * 309 + 1234 * (-1082) = 1
-----------------------------------
 4321 * (x - 309) + 1234 * (y - (-1082)) = 0 となるから

 4321 * (x - 309) = -1234 * (y + 1082) が成り立つ

 -1234 は x - 309 の約数で 4321 は y + 1082 の約数だから、t を整数とすると、

 x - 309 = -1234 * t => x = 309 - 1234 * t
 y + 1082 = 4321 * t => y = -1082 + 4321 * t

となる

それでは、実際に値を代入して (x, y) の組を求めてみましょう。

  t    x      y   z = 4321 * x + 1234 * y
 -----------------------------------------
 -4  5245 -18366  1
 -3  4011 -14045  1
 -2  2777  -9724  1
 -1  1543  -5403  1
  0   309  -1082  1
  1  -925   3239  1
  2 -2159   7560  1
  3 -3393  11881  1
  4 -4627  16202  1

プログラムも簡単です。Python でプログラムを作ると次のようになります。

リスト : 拡張ユークリッドの互除法

def ext_euclid(a, b):
    xs = a, 1, 0
    ys = b, 0, 1
    while ys[0] != 0:
        q, z = divmod(xs[0], ys[0])
        xs, ys = ys, (z, xs[1] - q * ys[1], xs[2] - q * ys[2])
    return xs

x, xa, xb をタプル xs に、y, ya, yb をタプル ys に、z, za, zb をタプルにまとめています。関数 divmod(x, y) は x と y の商と余りを求める関数です。あとはアルゴリズムをそのままプログラムしただけなので、難しいところはないと思います。それでは実際に試してみましょう。

>>> ext_euclid(8, 5)
(1, 2, -3)
>>> ext_euclid(9, 5)
(1, -1, 2)
>>> ext_euclid(11, 3)
(1, -1, 4)
>>> ext_euclid(4321, 1234)
(1, 309, -1082)
>>> ext_euclid(12357, 100102)
(1, 30127, -3719)

●水差し問題

拙作のページ Puzzle DE Programming に 水差し問題 がありますが、それを解くのに拡張ユークリッドの互除法を使うことができます。

たとえば、8 リットルの容器 A と 5 リットルの容器 B を使う場合、8 * x + 5 * y = gcd(8, 5) を解くと、8 * 2 - 5 * 3 = 1 という解を見つけることができます。これは A で 2 回汲み B で 3 回捨てると、1 リットルの水が残ることを表します。4 リットルの水を汲みだす場合、両辺を 4 倍して 8 * 8 - 5 * 12 = 4 になります。A で 8 回汲み B で 12 回捨ててもいいのですが、一般解を求めるともっと少ない回数の解を見つけることができます。

  8 * x + 5 *  y = 4
- 8 * 8 - 5 * 12 = 4
---------------------------------
  8 * (x - 8) + 5 * (y + 12) = 0

  x - 8 = -5 * t => x = 8 - 5 * t
  y + 12 = 8 * t => y = -12 + 8 * t

  t = 0, x = 8,  y = -12
  t = 1, x = 3,  y = -4
  t = 2, x = -2, y = 4

上記のように、B で 4 回汲み A で 2 回捨てる解を見つけることができます。

●合同式の徐算

拡張ユークリッドの互除法は、合同式の徐算で逆数を求めるときにも役に立ちます。a * b ≡ 1 (mod n) において a の逆数 b を求める場合、a * x + n * y = 1 を求めます。ここで、gcd(a, n) = 1 が条件になることに注意してください。両辺に mod n を施すと n * y は 0 になるので、a * x ≡ 1 (mod n) が成り立ちます。つまり、x が a の逆数になるわけです。簡単な例を示しましょう。

>>> for a in range(1, 11):
...     print(ext_euclid(a, 11))
...
(1, 1, 0)
(1, -5, 1)
(1, 4, -1)
(1, 3, -1)
(1, -2, 1)
(1, 2, -1)
(1, -3, 2)
(1, -4, 3)
(1, 5, -4)
(1, -1, 1)

n が素数の場合、n とその剰余は互いに素になります。合同式で x の負数 -x は (n - x) になるので、mod 11 における 1 から 10 までの逆数は次のようになります。

  x   1/x
 -------------------------
  0   ---
  1    1     1 mod 11 = 1
  2    6    12 mod 11 = 1
  3    4    12 mod 11 = 1
  4    3    12 mod 11 = 1
  5    9    45 mod 11 = 1
  6    2    12 mod 11 = 1
  7    8    56 mod 11 = 1
  8    7    56 mod 11 = 1
  9    5    45 mod 11 = 1
 10   10   100 mod 11 = 1

n が素数でない場合、n とその剰余は互いに素になるとは限りません。次の例を見てください。

>>> for x in range(1, 8): print(ext_euclid(x, 8))
...
(1, 1, 0)
(2, 1, 0)
(1, 3, -1)
(4, 1, 0)
(1, -3, 2)
(2, -1, 1)
(1, -1, 1)

2, 4, 6 と 8 の最大公約数は 1 にならないので、その逆数を求めることはできません。つまり、mod 8 の世界では 2, 4, 6 で割り算することはできないわけです。

また、フェルマーの小定理から逆数を求めることもできます。

[フェルマーの小定理]
n が素数で a と n が互いに素のとき an-1 ≡ 1 (mod n) が成り立つ

a * an-2 ≡1 (mod n) が成り立つので、逆数は an-2 になります。それでは実際に試してみましょう。

>>> for x in range(1, 11): print(x, x**9, (x**9) % 11)
...
1 1 1
2 512 6
3 19683 4
4 262144 3
5 1953125 9
6 10077696 2
7 40353607 8
8 134217728 7
9 387420489 5
10 1000000000 10

拡張ユークリッドの互除法と同じ結果になりましたね。なお、簡単な例ということで累乗を計算してから mod をとりましたが、実用的なプログラムではやらないほうがよいでしょう。mod n の世界で累乗を計算する関数を作ったほうがよいと思います。興味のある方は挑戦してみてください。

●参考 URL


2018 年 5 月

5月5日

●Python の高速化

Python のお話です。Python は使いやすいプログラミング言語ですが、実行速度はお世辞にも速いとは言えません。時間がかかる処理をC言語で記述してライブラリを作成し、それをインポートして呼び出す方法もありますが、Python のプログラムを変更なして高速化できるならば、その方が簡単で便利です。実際、そのようなライブラリに Numba があり、Numba をインポートすると Python に JIT (just in time) コンパイラを導入することができます。

Numba のインストールは pip を使うと簡単です。ところが、M.Hiroi の環境 (Windows10, Python 3.6.4) では、Numba のインストールには成功するのですが、Numba のインポートでエラーが発生します。対処方法がわからなかったので、今回は VirtualBox の Xubunts で Numba を試してみることにしました。Windows で Numba を使う場合、Anaconda のような科学技術計算向けのパッケージを導入したほうが簡単かもしれません。

pip は Python で書かれたパッケージ管理ツールです。次のコマンドで pip がインストールされているか確認することができます。

$ python3 -m pip -V
/usr/bin/python3: No module named pip

pip が入っていない場合は、次のコマンドでインストールしてください。

sudo apt install python3-pip

あとは pip で Numba をインストールするだけです。

$ python3 -m pip install numba

Numba の基本的な使い方は簡単で、高速化したい関数の前に @numba.jit を付けるだけです。これでその関数が JIT コンパイルされます。それでは実際に、いつもの「たらいまわし関数」で実行時間を計測してみましょう。

リスト : たらいまわし関数

import time, numba
                                                                                
def tak(x, y, z):
    if x <= y: return z
    return tak(tak(x - 1, y, z), tak(y - 1, z, x), tak(z - 1, x, y))

@numba.jit
def tak1(x, y, z):
    if x <= y: return z
    return tak1(tak1(x - 1, y, z), tak1(y - 1, z, x), tak1(z - 1, x, y))

s = time.time()
print(tak(22, 11, 0))
print(time.time() - s)
s = time.time()
print(tak1(22, 11, 0))
print(time.time() - s)
mhiroi@mhiroi-VirtualBox:~/python$ python3 tarai.py 
11
88.36137557029724
11
2.570927858352661

実行環境 : xubuntu 16.10 on VirtualBox, Intel Core i5-6200U 2.30GHz

素の Python3 では tak(22, 11, 0) に 88 秒もかかりますが、JIT コンパイルすることで 2.6 秒に短縮することができました。JIT コンパイラを使ったスクリプト言語には Node.js (JavaScript) や Julia などがありますが、それらの言語でたらいまわし関数を実行すると、Node.js で約 3 秒、Julia で 2 秒ほどかかります。これらの言語と同程度の速度を叩き出すのですから、Numba の JIT コンパイラは優秀ですね。Python のプログラムを高速化するときには試してみたいツールだと思ました。


2018 年 1 月

1月21日

●.Net Core

.NET のお話です。MicroSoft 社が開発したプログラミング言語 C# と F# は、.NET Framework の共通中間言語 (CLI) にコンパイルされて、共通言語ランタイム (CLR) 上で実行されます。.NET Framework 互換の環境にはオープンソースプロジェクトの「Mono」がありますが、このほかに、.NET Framework のサブセットとして MicroSoft 社がオープンソースで開発している .NET Core があります。.NET Core は Mono と同様にクロスプラットフォーム (Windows, Linux, MacOS) で動作します。

.NET Core の場合、.NET Core SDK をインストールすれば最低限の開発環境が整うので、C# や F# の学習が目的ならば .NET Core のほうが手軽で便利かもしれません。とりあえず Windows で試してみました。.NET Core SDK は .NET.NET Downloads からダウンロードすることができます。Windows の場合、インストーラをダウンロードして実行するだけです。

.NET Core SDK でプログラムを作る場合、dotnet というコマンドラインツール (CLI) を使います。たとえば、コンソールアプリケーションを作成する場合、プログラムを作成するディレクトリで次のコマンドを実行します。

dotnet new console -lang F#

dotnet new はプロジェクトを初期化するコマンドです。-lang は作成するプログラムの種類 (C# / F# / VB) を指定します。省略した場合は C# になります。これでカレントディレクトリにプログラム Program.fs とプロジェクトに必要なデータが生成されます。あとは Program.fs を書き換えて dotnet run と入力すれば、プログラムを実行することができます。

簡単といえば簡単なのですが、学習目的の小さなプログラムでプロジェクトを作るのは、ちょっと大げさなような気がします。まあ、今どきのツールはたいていそうなっているので、M.Hiroi の感覚がずれている (古い) のでしょう。F# にはインタプリタ (fsi.exe) があり、それを使えるとよかったのですが、起動方法がよくわかりませんでした。もしかしたら、.NET Core SDK では対応していないのかもしれません。


1月1日

あけましておめでとうございます

旧年中は大変お世話になりました
本年も M.Hiroi's Home Page をよろしくお願い申し上げます


Copyright (C) 2018 Makoto Hiroi
All rights reserved.

[ Home ]