M.Hiroi's Home Page

Memorandum

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

2018 年 12 月

12月1日

●格子点と直積集合

平面や空間などの座標において、各成分がすべて整数であるような点を「格子点 (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|...], [...|...]|...].

参考 URL

  1. 直積集合 - Wikepedia

2018 年 11 月

11月21日

●ラマヌジャンのタクシー数

笹川さんの 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) でも同様です。興味のある方は、いろいろなプログラミング言語で挑戦してみてください。

●追記 (2018/11/25)

最近のパソコンは高性能なので、単純な四重ループでも簡単に解くことができます。たとえば、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 倍ちょっと高速化することができました。なお、この結果にも重複解が含まれていることに注意してください。

●参考 URL

  1. 天才数学者ラマヌジャンのタクシー数の研究, (猫野さん)
  2. タクシー数 - Wikipedia

2018 年 10 月

10月2日

●Yahoo! ジオシティーズの閉鎖

このたび、Yahoo! ジオシティーズが 2019 年 3 月 31 日にサービスを終了することになりました。突然のことで大変驚いております。M.Hiroi's Home Page を開設したのは 2000 年で、かれこれ 18 年の長い間お世話になりました。閉鎖されるのは残念ですが、しかたがないですね。今まで本当にありがとうございました。

さて、M.Hiroi's Home Page の今後ですが、引っ越しして更新を続けていくつもりです。今回の件で転送サービスが用意されているのには助かりました。移転先が決まったらトップページでお知らせしなければと思ったのですが、しなくても大丈夫そうですね。それでは、これからも M.Hiroi's Home Page をよろしくお願い申し上げます。


2018 年 9 月

9月22日

●Easy-ISLisp の型推論

このたび、笹川さんが開発されている 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 を開発している笹川さんに感謝です。興味のある方は実際に試してみてください。


2018 年 6 月

6月24日

●有限体

●拡大体

Algorithms with Python: 拡張ユークリッドの互除法 へ移動


6月16日

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

Algorithms with Python: 拡張ユークリッドの互除法 へ移動


2018 年 5 月

5月5日

●Python の高速化

Python3 超入門: Python の高速化 へ移動


2018 年 1 月

1月21日

●.Net Core

.NET の説明は拙作のページ C# Programming, F# Programming をお読みください。


1月1日

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

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


Copyright (C) 2018 Makoto Hiroi
All rights reserved.

[ Home ]