平面や空間などの座標において、各成分がすべて整数であるような点を「格子点 (lattice point)」といいます。たとえば、二次元の座標 (x, y) で x と y の範囲が有限であれば、格子点は「直積集合 (direct product)」で求めることができます。直積は「デカルト積 (Cartesian product)」と呼ばれることもあります。
最近では、多くの処理系で集合演算をサポートしています。直積集合であれば、Gauche の cartesian-product、Julia の Iterators.product()、Python の itertools.product()、Ruby の product() などなど探せばいろいろあると思いますが、私たちでも簡単に作ることができます。
Python のように内包表記が使える処理系はとても簡単です。
リスト : 直積集合 (Python)
def product(xs, ys):
return [(x, y) for x in xs for y in ys]
end
>>> product([1, 2, 3], [4, 5, 6]) [(1, 4), (1, 5), (1, 6), (2, 4), (2, 5), (2, 6), (3, 4), (3, 5), (3, 6)]
三次元の座標は product() を 2 回使えばできるはずです。結果は次のようになりました。
>>> product(product([1, 2], [3, 4]), [5, 6]) [((1, 3), 5), ((1, 3), 6), ((1, 4), 5), ((1, 4), 6), ((2, 3), 5), ((2, 3), 6), ((2, 4), 5), ((2, 4), 6)] >>> product([1, 2], product1([3, 4], [5, 6])) [(1, (3, 5)), (1, (3, 6)), (1, (4, 5)), (1, (4, 6)), (2, (3, 5)), (2, (3, 6)), (2, (4, 5)), (2, (4, 6))]
集合 A, B の直積集合を A * B で表すと、一般に A * B = B * A は成り立ちません。また、(A * B) * C の要素は ((a, b), c) に、A * (B * C) の要素は (a, (b, c)) に、A * B * C の要素は (a, b, c) になるので、厳密に言えばこれらは異なる集合になります。
実際には、((a, b), c) や (a, (b, c)) を (a, b, c) と同一視して、複数の集合の直積を 2 つの集合の直積の繰り返しで定義してもよいようです。詳しくは 参考 URL 1 をお読みください。
これを Python でプログラムすると次のようになります。
リスト : 直積集合 (その二)
def product(*args):
if len(args) == 0: return [()]
ys = [(x,) for x in args[0]]
for i in range(1, len(args)):
ys = [(*y, x) for y in ys for x in args[i]]
return ys
最初に、先頭のリスト args[0] の要素をタプルに包み、それを格納したリストに変換します。あとは順番に、タプル y の後ろにリストの要素 x を追加していくだけです。Python の場合、リストやタプルの前に * を付けると、要素が展開されることに注意してください。
それでは実際に試してみましょう。
>>> product([1, 2, 3], [4, 5, 6]) [(1, 4), (1, 5), (1, 6), (2, 4), (2, 5), (2, 6), (3, 4), (3, 5), (3, 6)] >>> product([1, 2], [3, 4], [5, 6]) [(1, 3, 5), (1, 3, 6), (1, 4, 5), (1, 4, 6), (2, 3, 5), (2, 3, 6), (2, 4, 5), (2, 4, 6)] >>> product([1, 2], [3, 4], [5, 6], [7, 8]) [(1, 3, 5, 7), (1, 3, 5, 8), (1, 3, 6, 7), (1, 3, 6, 8), (1, 4, 5, 7), (1, 4, 5, 8), (1, 4, 6, 7), (1, 4, 6, 8), (2, 3, 5, 7), (2, 3, 5, 8), (2, 3, 6, 7), (2, 3, 6, 8), (2, 4, 5, 7), (2, 4, 5, 8), (2, 4, 6, 7), (2, 4, 6, 8)]
ご参考までに、Scheme (Gauche) と Prolog のプログラムを示します。どちらも product() の引数にはリストのリストを渡すことに注意してください。
リスト : 直積集合 (Scheme)
(define (flatmap f xs) (apply append (map f xs)))
(define (product xs)
(cond ((null? xs) (list '()))
((null? (cdr xs))
(map (lambda (x) (list x)) (car xs)))
(else
(flatmap (lambda (x) (map (lambda (ys) (cons x ys))
(product (cdr xs))))
(car xs)))))
gosh> (product '()) (()) gosh> (product '((1 2 3 4))) ((1) (2) (3) (4)) gosh> (product '((1 2 3) (4 5 6))) ((1 4) (1 5) (1 6) (2 4) (2 5) (2 6) (3 4) (3 5) (3 6)) gosh> (product '((1 2) (3 4) (5 6))) ((1 3 5) (1 3 6) (1 4 5) (1 4 6) (2 3 5) (2 3 6) (2 4 5) (2 4 6)) gosh> (product '((1 2) (3 4) (5 6) (7 8))) ((1 3 5 7) (1 3 5 8) (1 3 6 7) (1 3 6 8) (1 4 5 7) (1 4 5 8) (1 4 6 7) (1 4 6 8) (2 3 5 7) (2 3 5 8) (2 3 6 7) (2 3 6 8) (2 4 5 7) (2 4 5 8) (2 4 6 7) (2 4 6 8))
リスト : 直積集合 (Prolog) product([Xs], [X]) :- member(X, Xs). product([Xs | Ys], [X | Zs]) :- member(X, Xs), product(Ys, Zs).
?- product([], A). false. ?- product([[1,2,3]], A), write(A), fail. [1][2][3] false. ?- product([[1,2,3],[4,5,6]], A), write(A), fail. [1,4][1,5][1,6][2,4][2,5][2,6][3,4][3,5][3,6] false. ?- product([[1,2],[3,4],[5,6]], A), write(A), fail. [1,3,5][1,3,6][1,4,5][1,4,6][2,3,5][2,3,6][2,4,5][2,4,6] false. ?- product([[1,2],[3,4],[5,6],[7,8]], A), write(A), fail. [1,3,5,7][1,3,5,8][1,3,6,7][1,3,6,8][1,4,5,7][1,4,5,8][1,4,6,7][1,4,6,8] [2,3,5,7][2,3,5,8][2,3,6,7][2,3,6,8][2,4,5,7][2,4,5,8][2,4,6,7][2,4,6,8] false. ?- findall(A, product([[1,2],[3,4],[5,6],[7,8]], A), B). B = [[1, 3, 5, 7], [1, 3, 5, 8], [1, 3, 6, 7], [1, 3, 6, 8], [1, 4, 5, 7], [1, 4, 5|...], [1, 4|...], [1|...], [...|...]|...].
笹川さんの twitter に触発されて、M.Hiroi もタクシー数を求めるプログラムを書いてみました。ここでは問題を簡単にして、次式を満たす異なる整数 a, b, c, d を列挙することにします。
a ^ 3 + b ^ 3 = c ^ 3 + d ^ 3 (= x)
x が昇順になるように出力する場合、SWI-Prolog の CLPFD を使うと簡単です。
?- use_module(library(clpfd)). true. ?- A^3 + B^3 #= X, C^3 + D^3 #= X, | A #< B, C #< D, A #< C, D #< B, | [A, B, C, D] ins 1 .. 50, | labeling([min(X)], [A, B, C, D]), | write([X, A, B, C, D]), | nl, | fail. [1729,1,12,9,10] [4104,2,16,9,15] [13832,2,24,18,20] [20683,10,27,19,24] [32832,4,32,18,30] [39312,2,34,15,33] [40033,9,34,16,33] [46683,3,36,27,30] [64232,17,39,26,36] [65728,12,40,31,33] [110656,4,48,36,40] [110808,6,48,27,45] false.
labeling で min(X) を指定することで、最初に X が最小となる解を求めることができます。逆に、max(X) を指定すると最大値の解を求めることができます。このように簡単にプログラムできるのですが、SWI-Prolog の CLPFD は遅いので、ちょっと時間がかかるのが欠点です。他の処理系では、高速に解を求めることができるかもしれません。興味のある方は試してみてください。
ところで、整数 a, b, c, d の範囲を 1 から N として、条件を満たすものを列挙するだけでよければ、プログラムはぐっと簡単になります。たとえば、Prolog では次のようになります。
リスト : タクシー数
product(N, [X, A, B]) :-
between(1, N, A),
A1 is A + 1,
between(A1, N, B),
X is A ^ 3 + B ^ 3.
solver(N) :-
assert(taxi(0 ,0, 0)),
retractall(taxi(X, A, B)),
product(N, [X, A, B]),
assert(taxi(X, A, B)),
findall([C, D], taxi(X, C, D), Zs),
length(Zs, 2),
write([X, Zs]),
nl,
fail.
product() で整数の組を生成し、assert() でそれと三乗和を記憶します。そして、findall() で三乗和 X と等しい整数の組を求め、それが 2 つあれば write() で出力します。とても簡単ですね。
それでは実際に試してみましょう。
?- solver(50). [1729,[[1,12],[9,10]]] [4104,[[2,16],[9,15]]] [39312,[[2,34],[15,33]]] [40033,[[9,34],[16,33]]] [13832,[[2,24],[18,20]]] [32832,[[4,32],[18,30]]] [20683,[[10,27],[19,24]]] [64232,[[17,39],[26,36]]] [46683,[[3,36],[27,30]]] [110808,[[6,48],[27,45]]] [65728,[[12,40],[31,33]]] [110656,[[4,48],[36,40]]] false.
一瞬で 12 通りの解を求めることができました。solver(100) でも同様です。興味のある方は、いろいろなプログラミング言語で挑戦してみてください。
最近のパソコンは高性能なので、単純な四重ループでも簡単に解くことができます。たとえば、Julia でプログラムすると次のようになります。
リスト : タクシー数 (Julia)
# 単純な四重ループ版
function taxi4(n)
xs = Vector{Int}[]
for a = 1 : n, c = a + 1 : n, d = c + 1 : n, b = d + 1 : n
if a^3 + b^3 == c^3 + d^3
push!(xs, [a^3+b^3, a, b, c, d])
end
end
sort(xs, by=ys -> ys[1])
end
julia> taxi4(50)
12-element Array{Array{Int64,1},1}:
[1729, 1, 12, 9, 10]
[4104, 2, 16, 9, 15]
[13832, 2, 24, 18, 20]
[20683, 10, 27, 19, 24]
[32832, 4, 32, 18, 30]
[39312, 2, 34, 15, 33]
[40033, 9, 34, 16, 33]
[46683, 3, 36, 27, 30]
[64232, 17, 39, 26, 36]
[65728, 12, 40, 31, 33]
[110656, 4, 48, 36, 40]
[110808, 6, 48, 27, 45]
julia> @time taxi4(500)
2.524836 seconds (1.15 k allocations: 103.172 KiB)
566-element Array{Array{Int64,1},1}:
[1729, 1, 12, 9, 10]
[4104, 2, 16, 9, 15]
[13832, 2, 24, 18, 20]
[20683, 10, 27, 19, 24]
[32832, 4, 32, 18, 30]
[39312, 2, 34, 15, 33]
[40033, 9, 34, 16, 33]
[46683, 3, 36, 27, 30]
[64232, 17, 39, 26, 36]
[65728, 12, 40, 31, 33]
[110656, 4, 48, 36, 40]
[110808, 6, 48, 27, 45]
[134379, 12, 51, 38, 43]
⋮
[125461440, 238, 482, 336, 444]
[126217000, 250, 480, 295, 465]
[126516061, 165, 496, 262, 477]
[126906624, 124, 500, 307, 461]
[127062000, 235, 485, 330, 450]
[133139727, 228, 495, 306, 471]
[137854899, 243, 498, 395, 424]
[138090771, 294, 483, 392, 427]
[143604279, 359, 460, 408, 423]
[143992576, 292, 492, 414, 418]
[145908847, 346, 471, 360, 463]
[147748104, 306, 492, 384, 450]
さすがに単純な四重ループでは、整数の範囲を増やすと遅くなります。なお、この結果には重複解が含まれていることに注意してください。
[87539319, 167, 436, 228, 423] [87539319, 167, 436, 255, 414] [87539319, 228, 423, 255, 414] [119824488, 11, 493, 90, 492] [119824488, 11, 493, 346, 428] [119824488, 90, 492, 346, 428]
これらは条件 a^3 + b^3 = c^3 + d^3 = e^3 + f^3 を満たしています。この条件を満たす最小の数が 87,539,319 です。
次に示すように辞書を使って二重ループにすると速くなります。
リスト : タクシー数 (2)
# 二重ループ版
function taxi2(n)
xs = []
table = Dict{Int, Tuple{Int, Int}}()
for a = 1 : n, b = a + 1 : n
x = a^3 + b^3
if haskey(table, x)
push!(xs, [x, (a, b), table[x]])
else
table[x] = (a, b)
end
end
sort(xs, by=ys -> ys[1])
end
julia> @time taxi2(500)
0.221643 seconds (368.06 k allocations: 27.129 MiB, 3.21% gc time)
564-element Array{Any,1}:
Any[1729, (9, 10), (1, 12)]
Any[4104, (9, 15), (2, 16)]
Any[13832, (18, 20), (2, 24)]
Any[20683, (19, 24), (10, 27)]
Any[32832, (18, 30), (4, 32)]
Any[39312, (15, 33), (2, 34)]
Any[40033, (16, 33), (9, 34)]
Any[46683, (27, 30), (3, 36)]
Any[64232, (26, 36), (17, 39)]
Any[65728, (31, 33), (12, 40)]
Any[110656, (36, 40), (4, 48)]
Any[110808, (27, 45), (6, 48)]
Any[134379, (38, 43), (12, 51)]
⋮
Any[125461440, (336, 444), (238, 482)]
Any[126217000, (295, 465), (250, 480)]
Any[126516061, (262, 477), (165, 496)]
Any[126906624, (307, 461), (124, 500)]
Any[127062000, (330, 450), (235, 485)]
Any[133139727, (306, 471), (228, 495)]
Any[137854899, (395, 424), (243, 498)]
Any[138090771, (392, 427), (294, 483)]
Any[143604279, (408, 423), (359, 460)]
Any[143992576, (414, 418), (292, 492)]
Any[145908847, (360, 463), (346, 471)]
Any[147748104, (384, 450), (306, 492)]
10 倍ちょっと高速化することができました。なお、この結果にも重複解が含まれていることに注意してください。
このたび、Yahoo! ジオシティーズが 2019 年 3 月 31 日にサービスを終了することになりました。突然のことで大変驚いております。M.Hiroi's Home Page を開設したのは 2000 年で、かれこれ 18 年の長い間お世話になりました。閉鎖されるのは残念ですが、しかたがないですね。今まで本当にありがとうございました。
さて、M.Hiroi's Home Page の今後ですが、引っ越しして更新を続けていくつもりです。今回の件で転送サービスが用意されているのには助かりました。移転先が決まったらトップページでお知らせしなければと思ったのですが、しなくても大丈夫そうですね。それでは、これからも M.Hiroi's Home Page をよろしくお願い申し上げます。
このたび、笹川さんが開発されている Easy-ISLisp (EISL) の最新版 (ver 0.90) に、型推論とデータ型を指定する機能が追加されましたので、VirtualBox の xubuntu 上で試してみました。EISL は下記のページからダウンロードすることができます。
コンパイラの使い方や型推論については下記のページを参照してください。
EISL (Linux 版) は TAR でまとめられているので、tar -xvf で展開してください。その中のファイル compile.o がコンパイラのオブジェクトファイルです。 ./eisl -c compile.o でオブジェクトファイルをリンクしてください。これでコンパイラの動作が高速になります。
それではいつもの「たらいまわし関数」をコンパイルしてみましょう。
リスト : たらいまわし関数 (tak.lsp)
(defun tak (x y z)
(if (<= x y)
z
(tak (tak (- x 1) y z)
(tak (- y 1) z x)
(tak (- z 1) x y))))
コンパイルは簡単で (compile-file "tak.lsp") とするだけです。これでオブジェクトファイル tak.o が生成されます。動的リンクも簡単です。(load "tak.o") とするだけです。実行結果は次のようになりました。
$ ./eisl -c compiler.o Easy-ISLisp Ver0.90 > (compile-file "tak.lsp") type inference initialize pass1 pass2 compiling TAK finalize invoke GCC T > (load "tak.o") T > (time (tak 22 11 0)) Elapsed Time(second)=17.031096 <undef>
前回試してみた ver 0.81 よりも少し遅くなっているようです。そこで、引数のデータ型を小整数 (fixnum) に指定して試してみます。プログラムは次のようになります。
リスト : たらいまわし関数 (tak1.lsp)
(defun tak (x y z)
(the <fixnum> x) (the <fixnum> y) (the <fixnum> z)
(if (<= x y)
z
(tak (tak (- x 1) y z)
(tak (- y 1) z x)
(tak (- z 1) x y))))
EISL の場合、データ型の指定には the という構文を使います。これで引数 x, y, z のデータ型が fixnum であることをコンパイラに教えることができます。実行結果は次のようになりました。
$ ./eisl -c compiler.o Easy-ISLisp Ver0.90 > (compile-file "tak1.lsp") type inference initialize pass1 pass2 compiling TAK finalize invoke GCC T > (load "tak1.o") T > (time (tak 22 11 0)) Elapsed Time(second)=1.329383 <undef>
10 倍以上高速化することができました。SBCL で最適化を施すと約 1.67 秒になるので、それよりも少し速い結果になりました。ここまで速くなるとは M.Hiroi も大変驚きました。EISL を開発している笹川さんに感謝です。興味のある方は実際に試してみてください。
数学の世界では、加減乗除ができる数の集合のことを「体」といいます。特に、数が 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 次方程式を列挙すると次のようになります。
この中で解を持たないのは式 4 です。『x の関数 f(x) において、f(a) = 0 ならば (x - a) で割り切れる』という性質を「因数定理」といいます。式 1, 2, 3 を因数分解すると次のようになります。
式 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)
変数の個数が式の本数よりも多い方程式を「不定方程式」といいます。一次不定方程式 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 の世界で累乗を計算する関数を作ったほうがよいと思います。興味のある方は挑戦してみてください。
Python のお話です。Python は使いやすいプログラミング言語ですが、実行速度はお世辞にも速いとは言えません。時間がかかる処理をＣ言語で記述してライブラリを作成し、それをインポートして呼び出す方法もありますが、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 のプログラムを高速化するときには試してみたいツールだと思ました。
.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 では対応していないのかもしれません。