分数を小数に直すことを考えます。1 / 2, 1 / 5, 1 / 8, 1 / 10 などは 0.5, 0.2, 0.125, 0.1 のように途中で割り切れて、有限な桁で表すことができます。これを「有限小数」といいます。ところが、1 / 3, 1 / 6, 1 / 7 は、次のように有限な桁では表すことができません。
1 / 3 = 0.333333 .., 1 / 6 = 0.166666 ... 1 / 7 = 0.142857142857142857 ...
1 / 3 は 3 を無限に繰り返し、1 / 6 は 0.1 のあと 6 を無限に繰り返し、1 / 7 は数列 142857 を無限に繰り返します。このような少数を「循環小数 (repeating decimal) 」といい、繰り返される数字の列を「循環節」といいます。分数を小数に直すと、有限小数か循環小数のどちらかになります。
循環小数は次のように循環節の始めと終わりを点で示します。
.
1 / 3 = 0.3
.
1 / 6 = 0.16
. .
1 / 7 = 0.142857
このほかに、循環節に下線を引いたり、カッコで囲む表記法もあります。
分数を循環小数に直すのは簡単です。筆算と同じように計算していくだけです。次の図を見てください。
0 1 4 2 8 5 7
----------------
7 ) 1
0
-----
1 0 <-- 余り 1
7
-------
3 0
2 8
-------
2 0
1 4
-------
6 0
5 6
-------
4 0
3 5
-------
5 0
4 9
-----
1 <-- 余り 1
7 で割った余り 1 が 2 回現れていますね。これから先は同じ数列を繰り返すことがわかります。7 の剰余は 0 から 6 まであり、剰余が 0 の場合は割り切れるので循環小数にはなりません。現れる余りの数が限られているので、割り切れなければ循環することになるわけです。また、このことから循環節の長さは分母の値よりも小さいことがわかります。
それでは Haskell でプログラムを作ってみましょう。
リスト : 循環小数を求める
repdec :: Integer -> Integer -> ([Integer], [Integer])
repdec m n = iter m [] []
where iter m xs ys =
if q == 0
then (reverse (p:ys), [0])
else if q `elem` xs
then splitAt (length a + 1) (reverse (p:ys))
else iter (q * 10) (q:xs) (p:ys)
where p = m `div` n -- 商
q = m `mod` n -- 余
a = takeWhile (/= q) (reverse xs)
関数 repdec m n は m / n を循環小数に変換します。返り値は 2 つのリストを格納したタプルで、先頭のリストが冒頭の小数部分を、次のリストが循環節の部分を表します。途中で割り切れる場合は循環節を [0] とします。これ以降、冒頭の小数部分を有限小数部分と記述することにします。
実際の処理は局所関数 iter で行います。第 2 引数 xs が余りを格納するリスト、第 3 引数 ys が商を格納するリストです。最初に m と n の商と余りを計算して、変数 p と q にセットします。q が 0 ならば割り切れたので有限小数です。p:ys を reverse で反転し、それと [0] をタプルに格納して返します。
割り切れない場合、余り q が xs にあるか elem でチェックします。同じ値が無ければ iter を再帰呼び出しします。同じ値を見つけた場合、reverse で xs を反転し、takeWhile で先頭から q の手前までの要素を取り出し、変数 a にセットします。有限小数部分は a の長さ +1 までなので、splitAt で商のリスト (reverse (p:ys)) を分割して返します。
それでは実行してみましょう。
*Main> repdec 1 3 ([0],[3]) *Main> repdec 1 6 ([0,1],[6]) *Main> repdec 1 7 ([0],[1,4,2,8,5,7]) *Main> repdec 1 2 ([0,5],[0]) *Main> repdec 1 9 ([0],[1]) *Main> repdec 1 11 ([0],[0,9]) *Main> repdec 355 113 ([3],[1,4,1,5,9,2,9,2,0,3,5,3,9,8,2,3,0,0,8,8,4,9,5,5,7,5,2,2,1,2,3,8,9,3,8,0,5, 3,0,9,7,3,4,5,1,3,2,7,4,3,3,6,2,8,3,1,8,5,8,4,0,7,0,7,9,6,4,6,0,1,7,6,9,9,1,1,5, 0,4,4,2,4,7,7,8,7,6,1,0,6,1,9,4,6,9,0,2,6,5,4,8,6,7,2,5,6,6,3,7,1,6,8])
正常に動作していますね。
循環小数を分数に直すことも簡単にできます。循環小数 - Wikipedia によると、有限小数部分を a、循環節の小数表記を b、節の長さを n とすると、循環小数 x は次の式で表すことができる、とのことです。
1 1 1
x = a + b * (1 + ------ + ------- + ------- + ...)
10^n 10^2n 10^3n
10^n
= a + b * ----------
10^n - 1
ここで、カッコの中は初項 1, 公比 1 / (10^n) の無限等比級数になるので、その和は次の公式より求めることができます。
∞ a Σ arn = ----- (ただし |r| < 1) n=0 1 - r
簡単な例を示しましょう。
. .
0.142857 =
10^6
0.142857 * -------- = 142857 / 999999 = 1 / 7
10^6 - 1
.
0.16 =
10 1 6 10 15
0.1 + 0.06 * ---- = -- + --- * -- = ----- = 1 / 6
9 10 100 9 90
プログラムを作る場合、次のように考えると簡単です。
有限小数部分を表すリストを xs とすると
有限小数部分 = p / q
ただし p = xs を整数に変換
q = 10 ^ (length xs - 1)
循環節を表すリストを ys とすると
循環節 = p' / q'
ただし p' = ys を整数に変換
q' = 10 ^ (length ys)
p p' p * q' + p'
循環小数 = --- + -------- = -------------
q q * q' q * q'
冒頭の有限小数部分を分数に変換するのは簡単ですね。先頭が整数部分になるので、小数部分の桁はリスト xs の長さ - 1 になります。循環節だけを分数に変換する場合、たとえば 1 / 7 の循環節は [1,4,2,8,5,7] になりますが、分子 p' は 0.142857 * 10^6 = 142857 となるので、循環節を表すリストを整数に変換するだけでよいことがわかります。有限小数部分がある場合は、その桁数だけ循環節の部分を右シフトすればよいので、p' / q' に 1 / q を掛けるだけです。
これを Haskell でプログラムすると、次のようになります。
リスト : 循環小数を分数に直す
fromRepdec :: ([Integer], [Integer]) -> (Integer, Integer)
fromRepdec (xs, ys) = (b `div` c, a `div` c)
where
toNum = foldl (\a x -> a * 10 + x) 0
-- 有限小数の部分を分数に直す
p = toNum xs
q = 10 ^ (fromIntegral (length xs - 1))
-- 循環節を分数に直す
p' = toNum ys
q' = (10 ^ (fromIntegral (length ys))) - 1
-- 有限小数 + 循環節
a = q * q'
b = q' * p + p'
c = gcd a b
アルゴリズムをそのままプログラムしただけなので、とくに難しいところは無いと思います。
それでは実行してみましょう。
*Main> fromRepdec ([0], [1]) (1,9) *Main> fromRepdec ([0], [3]) (1,3) *Main> fromRepdec ([0], [1,4,2,8,5,7]) (1,7) *Main> fromRepdec ([0,1], [6]) (1,6) *Main> fromRepdec ([0,1], [0]) (1,10) *Main> fromRepdec ([0,2,5], [0]) (1,4) *Main> fromRepdec ([0], [1,2]) (4,33) *Main> fromRepdec ([0], [1,2,3]) (41,333) *Main> fromRepdec ([0], [1,2,3,4]) (1234,9999) *Main> fromRepdec ([0], [1,2,3,4,5]) (4115,33333) *Main> fromRepdec ([0], [1,2,3,4,5,6]) (41152,333333) *Main> fromRepdec ([0], [1,2,3,4,5,6,7]) (1234567,9999999) *Main> repdec 1234567 9999999 ([0],[1,2,3,4,5,6,7])
1234567 / 9999999 の循環節が分子と同じ 1234567 になるのは面白いですね。
整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。簡単な例を示します。
n = 6
6 分割 : 1 + 1 + 1 + 1 + 1 + 1
5 分割 : 1 + 1 + 1 + 1 + 2
4 分割 : 1 + 1 + 1 + 3
1 + 1 + 2 + 2
3 分割 : 1 + 1 + 4
1 + 2 + 3
2 + 2 + 2
2 分割 : 1 + 5
2 + 4
3 + 3
1 分割 : 6
6 の場合、分割の仕方は 11 通りあります。この数を「分割数」といいます。
分割数は拙作のページ Yet Another Haskell Problems 問題 57 で出題しました。解答 57 では、ナイーブな方法と動的計画法を使ったプログラムを示しましたが、数が大きくなると動的計画法を使ったプログラムでも遅くなります。実際に試してみましょう。
*Main> partition_number' 1000 24061467864032622473692149727991 (0.27 secs, 103563064 bytes) *Main> partition_number' 2000 4720819175619413888601432406799959512200344166 (1.23 secs, 429698732 bytes) *Main> partition_number' 4000 1024150064776551375119256307915896842122498030313150910234889093895 (5.37 secs, 1669094332 bytes) 実行環境 : Windows 7, Core i7-2670QM 2.20GHz, GHCi, version 7.4.1
分割数 - Wikipedia によると、次の漸化式を使うと分割数を高速に求めることができるそうです。
p(k) = p(k - 1) + p(k - 2) - p(k - 5) - p(k - 7) + p(k - 12) + p(k - 15) - p(k - 22) - ...
漸化式の説明を 分割数 - Wikipedia より引用します。
『ここで p(0) = 1 および負の整数 k に対して p(k) = 0 とし、和は (1/2)n(3n - 1) の形（ただし n は正または負の整数全体を走る）の一般五角数全体にわたってとるものとする（順に n = 1, -1, 2, -2, 3, -3, 4, -4 ..., とすると、値として 1, 2, 5, 7, 12, 15, 22, 26, 35, 40, 51, ... が得られる）。和における符号は交互に +, +, -, -, +, +, ... と続く。』
分割数 p(k) は k - 1 以下の分割数がわかれば求めることができます。この漸化式も動的計画法を使えば簡単にプログラムできます。次のリストを見てください。
リスト : 分割数
import qualified Data.Map as M
import Data.Maybe
-- 五角数
pentagon :: Integer -> Integer
pentagon x = x * (3 * x - 1) `div` 2
-- 分割数を求める
partNum :: M.Map Integer Integer -> Integer -> Integer -> Integer -> Integer
partNum m k n a =
if p1 >= 0 && p2 >= 0 then partNum m k (n + 2) (a + x1 + x2)
else if p1 >= 0 then a + x1
else a
where p1 = k - pentagon n
p2 = k - pentagon (-n)
x1 = fromJust $ M.lookup p1 m
x2 = fromJust $ M.lookup p2 m
partitionNum :: Integer -> Integer
partitionNum x = iter 2 $ M.fromList [(0,1), (1,1)]
where iter n m
| n > x = fromJust $ M.lookup x m
| otherwise = let v = partNum m n 1 0 - partNum m n 2 0
in iter (n + 1) (M.insert n v m)
分割数 p(k) を記憶するために Data.Map を使います。漸化式の計算は関数 partNum で行います。partNum の第 1 引数は Map で、第 2 引数 k の分割数を求めます。n は漸化式の項の番号を表し、第 4 引数の a が累積変数です。
整数 n = 1, -1, 3, -3, 5, -5 ... に対応する項の符号は + で、n = 2, -2, 4, -4, 6, -6 ... に対応する項の符号は - になります。引数 n を 1 から開始して n を +2 ずつ増やしていけば、+ の項の合計値を求めることができます。n を 2 から開始すれば - の項の合計値を求めることができます。分割数は + の項の値から - の項の値を引き算するだけです。
あとは、局所関数 iter で n を 2 から x まで順番に増やして n の分割数を求め、それをマップ m に格納していくだけです。最後に x の値を返します。
それでは実行してみましょう。
*Main> partition_number' 5000 169820168825442121851975101689306431361757683049829233322203824652329144349 (8.91 secs, 2629616996 bytes) *Main> partitionNum 5000 169820168825442121851975101689306431361757683049829233322203824652329144349 (1.45 secs, 185374304 bytes) *Main> partitionNum 10000 36167251325636293988820471890953695495016030339315650422081868605887952568754066 420592310556052906916435144 (4.07 secs, 530398016 bytes) *Main> partitionNum 50000 36261860971416678445921408915956337281653830825277850490158727554141099042567120 82718122747316610565824630881772910217544261659239432670671532413858378256188987 33387712189158660795738975053844747471259297926371901246185871979162730248973954 8263 (48.20 secs, 6412401312 bytes) 実行環境 : Windows 7, Core i7-2670QM 2.20GHz, GHCi, version 7.4.1
配列を使うと速くなるだろうと思って試してみたところ、マップよりも少しだけ速くなりました。コンパイルするとその差はもっと大きくなるかもしれません。興味のある方はいろいろ試してみてください。
ご参考までに配列を使ったプログラムを示します。
リスト : 分割数 (配列版)
import Data.Array.IO
-- 分割数を求める
partNum' :: IOArray Integer Integer -> Integer -> Integer -> Integer -> IO Integer
partNum' m k n a = do
let p1 = k - pentagon n
p2 = k - pentagon (-n)
if p1 < 0 then return a
else do
x1 <- readArray m p1
if p2 < 0
then return $ a + x1
else do x2 <- readArray m p2
partNum' m k (n + 2) (a + x1 + x2)
partitionNum' :: Integer -> IO Integer
partitionNum' x = do
a <- newArray_ (0, x) :: IO (IOArray Integer Integer)
writeArray a 0 1
writeArray a 1 1
iter 2 a
where iter n a
| n > x = readArray a x
| otherwise = do v1 <- partNum' a n 1 0
v2 <- partNum' a n 2 0
let v = v1 - v2
v `seq` writeArray a n v
iter (n + 1) a
点を多角形の形に並べたとき、その総数を「多角数 (polygonal number) 」といいます。三角形に配置したものを三角数、四角形に配置したものを四角数、五角形に配置したものを五角数といいます。三角数、四角数、五角数を下図に示します。
１ ３ ６ 10 15
● ● ● ● ●
●● ●● ●● ●●
●●● ●●● ●●●
●●●● ●●●●
●●●●●
図 : 三角数
１ ４ ９ 16 25
● ●● ●●● ●●●● ●●●●●
●● ●●● ●●●● ●●●●●
●●● ●●●● ●●●●●
●●●● ●●●●●
●●●●●
図 : 四角数
１ ５ 12 22
● ● ● ●
● ● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ● ● ●
● ● ●
● ● ● ●
● ● ● ●
● ● ● ●
図 : 五角数
多角数の点の増分を表に示すと、次のようになります。
n 三角数 四角数 五角数
---+-----------------------------------------------------------
1 | 1 1 1
2 | 3 = 1+2 4 = 1+3 5 = 1+4
3 | 6 = 1+2+3 9 = 1+3+5 12 = 1+4+7
4 | 10 = 1+2+3+4 16 = 1+3+5+7 22 = 1+4+7+10
5 | 15 = 1+2+3+4+5 25 = 1+3+5+7+9 35 = 1+4+7+10+13
6 | 21 = 1+2+3+4+5+6 36 = 1+3+5+7+9+11 51 = 1+4+7+10+13+16
・・・・・・ ・・・・・・・ ・・・・・
n | n(n + 1) / 2 n^2 n(3n - 1) / 2
表を見ればお分かりのように、三角数は公差 1、四角数は公差 2、五角数は公差 3、p 角数は公差 p - 2 の等差数列の和になります。初項を a, 公差を d とすると、等差数列の和 Sn は次式で求めることができます。
Sn = n(2a + (n - 1)d) / 2
a = 1, d = p - 2 を代入して計算すると、多角数 Pp,n は次式で求めることができます。
Pp,n = ((p - 2)n^2 - (p - 4)n) / 2
この式を Haskell でプログラムすると、次のようになります。
リスト : 多角数 polygonalNum :: Integer -> Integer -> Integer polygonalNum p n = ((p - 2) * n * n - (p - 4) * n) `div` 2
それでは実行してみましょう。
*Main> map (polygonalNum 3) [1..20] [1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,136,153,171,190,210] *Main> map (polygonalNum 4) [1..20] [1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400] *Main> map (polygonalNum 5) [1..20] [1,5,12,22,35,51,70,92,117,145,176,210,247,287,330,376,425,477,532,590] *Main> map (polygonalNum 6) [1..20] [1,6,15,28,45,66,91,120,153,190,231,276,325,378,435,496,561,630,703,780] *Main> map (polygonalNum 7) [1..20] [1,7,18,34,55,81,112,148,189,235,286,342,403,469,540,616,697,783,874,970] *Main> map (polygonalNum 8) [1..20] [1,8,21,40,65,96,133,176,225,280,341,408,481,560,645,736,833,936,1045,1160]
p 角数の初項は 1 で、第 2 項は p になります。多角数 - Wikipedia によると、多角数にはいろいろな面白い性質があるようです。
整数 m が p 角数か判定するプログラムも簡単です。次式で n が整数値になれば、m は p 角数の第 n 項であることがわかります。
(p - 4) + √((p - 4)^2 + 8(p - 2)m)
n = -----------------------------------
2(p - 2)
リスト : 多角数の判定
isPolygonalNum :: Integer -> Integer -> Bool
isPolygonalNum p m =
let a = p - 2
b = p - 4
c = b * b + 8 * a * m
d = floor (sqrt (fromIntegral c))
in d * d == c && (b + d) `mod` (2 * a) == 0
p - 2 を変数 a に、p - 4 を変数 b に、ルートの中を計算して変数 c にセットします。sqrt で c の平方根を求め、floor で整数値に変換して d にセットします。あとは、d * d が c に等しく、b + d が 2 * a で割り切れることを確認するだけです。
それでは実行してみましょう。
*Main> isPolygonalNum 3 55 True *Main> isPolygonalNum 3 54 False *Main> isPolygonalNum 4 64 True *Main> isPolygonalNum 4 65 False *Main> isPolygonalNum 5 117 True *Main> isPolygonalNum 5 118 False
もちろん、n の値を求めることもできます。興味のある方はプログラムを修正してみてください。
N 個の要素の順列は N! 通りあります。この順列に 0 から N! - 1 までの番号をつけることを考えます。拙作のページ Puzzle DE Programming 幅優先探索の高速化 -- ８パズルを解く -- では、配列を使った方法を紹介しましたが、リストを使ってプログラムすることもできます。
たとえば、0 から 8 までの 9 個の整数の順列で、番号の振り方を考えてみましょう。最初が 0 で始まるパターンは 8! = 40320 通りありますね。このパターンには 0 - 40319 までの番号を割り当てます。そして、1 で始まるパターンには 40320 - 80639 までの番号を割り当てます。残りのパターンも同じです。
次に 2 番目の数字を考えましょう。01 で始まるパターンは 7! = 5040 通りあります。 したがって、このパターンには 0 - 5039 までの番号を割り当てます。10 で始まるパターンには 40320 - 45359 までの番号を、12 で始まるパターンには 45360 - 50399 までの番号を割り当てます。あとはこれを 9 番目までの数字まで続ければ、すべてパターンに番号を割り当てることができます。
では実際に 867254301 というパターンで試してみましょう。次の図を見てください。
８ 8 * 8! ６ [0 1 2 3 4 5 6 7] : 8*8! + 6*7! ７ [0 1 2 3 4 5 7] : 8*8! + 6*7! + 6*6! ２ [0 1 2 3 4 5] : 8*8! + 6*7! + 6*6! + 2*5! ５ [0 1 3 4 5] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! ４ [0 1 3 4] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! ３ [0 1 3] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! + 2*2! ０ [0 1] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! + 2*2! + 0*1! １ [1] : 番号：357478
注意すべき点は、数字をそのまま掛け算してはいけないところです。たとえば、7 に注目してください。このとき、残されている数字は 0, 1, 2, 3, 4, 5, 7 がありますね。番号は順番に振っていくので、867 は 86 で始まるパターンの 6*6! 番目から始まるのです。つまり、残っている数字の中で何番目に位置しているのかを求める必要があります。
それでは、Haskell でプログラムを作りましょう。次のリストを見てください。
リスト : 順列を番号に変換する
import Data.List
import Data.Maybe
fromPerm :: Eq a => [a] -> [a] -> Integer
fromPerm xs ys = iter xs ys 0
where
iter [_] _ n = n
iter zs@(x:xs) ys n = iter xs (delete x ys) (n + m * (fromIntegral i))
where m = fact $ fromIntegral (length zs - 1)
i = fromJust $ findIndex (== x) ys
関数 fromPerm の第 1 引数が番号に変換する順列、第 2 引数が要素を昇順に並べたリストです。実際の処理は局所関数 iter で行います。iter の第 3 引数は累積変数で、ここに順列の番号が求まります。
第 1 引数の要素が一つになったら順列を番号に変換できたので、第 3 引数の n を返します。そうでなければ、第 1 引数 zs の先頭要素 x に番号を割り当てます。最初に、zs の長さから 1 を引いた値の階乗を求めて、変数 m にセットします。次に、findIndex で ys にある x の位置を求めて、変数 i にセットします。findIndex の返り値は Maybe 型です。この場合、x は必ず見つかるので formJust で Just から値を取り出しています。あとは n に i * m を加算するだけです。iter を再帰呼び出しするときは、delete で ys から x を削除することをお忘れなく。
それでは実行してみましょう。
*Main> fromPerm [1,2,3,4] [1,2,3,4] 0 *Main> fromPerm [2,1,3,4] [1,2,3,4] 6 *Main> fromPerm [3,2,1,4] [1,2,3,4] 14 *Main> fromPerm [4,3,2,1] [1,2,3,4] 23 *Main> map (\xs -> fromPerm xs [1..4]) $ permutation 4 [1,2,3,4] [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] *Main> fromPerm [0,1,2,3,4,5,6,7,8] [0..8] 0 *Main> fromPerm [4,5,6,7,8,0,1,2,3] [0..8] 184896 *Main> fromPerm [8,7,6,5,4,3,2,1,0] [0..8] 362879
permutation は拙作のページ Haskell 入門 順列と組み合わせ で作成した順列を生成する関数です。fromPerm は正常に動作していますね。
番号を順列に変換することも簡単です。次のリストを見てください。
リスト : 番号を順列に変換する (番号は 0 から開始)
toPerm :: Eq a => Integer -> [a] -> [a]
toPerm n xs = iter n xs []
where
iter _ [] a = reverse a
iter n xs a =
iter (n - m * p) (delete x xs) (x : a)
where m = fact (fromIntegral(length xs) - 1)
p = n `div` m
x = xs !! (fromIntegral p)
toPerm の第 1 引数が順列を表す番号、第 2 引数が要素を昇順に並べたリストです。実際の処理は局所関数 iter で行います。iter の第 3 引数は累積変数で、ここに順列が求まります。
iter の第 2 引数が空リストになったら番号を順列に変換できたので、第 3 引数 a を reverse で反転して返します。そうでなければ、xs の長さから 1 を引いた値の階乗を求めて変数 m にセットします。n `div` m で xs にある要素の位置がわかるので、それを演算子 !! で取り出して、変数 x にセットします。あとは、iter を再帰呼び出しするとき、n から m * p を引き算し、delete で x を xs から取り除き、x を累積変数 a に追加します。これで番号に対応した順列を求めることができます。
それでは実行してみましょう。
*Main> toPerm 0 [1,2,3,4] [1,2,3,4] *Main> toPerm 12 [1,2,3,4] [3,1,2,4] *Main> toPerm 23 [1,2,3,4] [4,3,2,1] *Main> map (\x -> toPerm x [1,2,3,4]) [0..23] [[1,2,3,4],[1,2,4,3],[1,3,2,4],[1,3,4,2],[1,4,2,3],[1,4,3,2],[2,1,3,4],[2,1,4,3], [2,3,1,4],[2,3,4,1],[2,4,1,3],[2,4,3,1],[3,1,2,4],[3,1,4,2],[3,2,1,4],[3,2,4,1], [3,4,1,2],[3,4,2,1],[4,1,2,3],[4,1,3,2],[4,2,1,3],[4,2,3,1],[4,3,1,2],[4,3,2,1]] *Main> toPerm 0 [0..8] [0,1,2,3,4,5,6,7,8] *Main> toPerm 362879 [0..8] [8,7,6,5,4,3,2,1,0]
正常に動作していますね。
ピタゴラスの定理は、平面幾何学において直角三角形の辺を a, b, c (a + b > c) とすると、次式が成り立つという皆さんお馴染みの有名な定理です。
a^2 + b^2 = c^2
ピタゴラス数 (Pythagoras triple) は、上式を満たす自然数の組 (a, b, c) のことで、とくに a, b, c が互いに素であるとき、(a, b, c) を原始ピタゴラス数といいます。たとえば、(3, 4, 5), (5, 12, 13), (8, 15, 17) などがあります。
ピタゴラスの定理 - Wikipedia によると、原始ピタゴラス数は次の方法で簡単に見つけることができるそうです。
本ページでは証明を割愛しますが、inamori さんのブログ 桃の天然水 ピタゴラス数 の説明がわかりやくて参考になると思います。inamori さんに感謝いたします。
それでは、実際にプログラムを作ってみましょう。使用するプログラミング言語は Haskell です。三辺の合計値 (a + b + c) が k 以下の原始ピタゴラス数をすべて求めることにします。プログラムは次のようになります。
リスト : 三辺の合計値が k 以下の原始ピタゴラス数を求める
primitivePythagoras :: Integer -> [(Integer, Integer, Integer)]
primitivePythagoras k =
foldr (\m xs -> iter 1 m ++ xs) [] [2 .. mh]
where mh = ceiling (sqrt (fromIntegral k * 0.5))
iter n m
| n >= m || 2 * m * (m + n) > k = []
| odd (m - n) && gcd m n == 1 =
let a = m * m - n * n
b = 2 * m * n
c = m * m + n * n
in (if a < b then (a, b, c) else (b, a, c)) : iter (n + 1) m
| otherwise = iter (n + 1) m
三辺の合計値は次の式で求めることができます。
a + b + c = 2mn + m^2 - n^2 + m^2 + n^2
= 2m(m + n)
上式から三辺の合計値は必ず偶数になることがわかります。変数 m の値を √(k/2) とすると、n = 1 のときに三辺の合計値は 2m^2 + 2m = k + 2 * √(k/2) となり、k よりも大きくなるので、m の上限値を √(k/2) とします。この値を変数 mh にセットします。ceiling で小数点を切り上げていることに注意してください。
foldr で m の値を 2 から mh まで増やしていき、ラムダ式の中の局所関数 iter で原始ピタゴラス数を生成します。n が m 以上になるか、三辺の合計値が k を超えたならば空リストを返します。そうでなければ、m - n が奇数で、m と n が互いに素であることを確認します。そして、辺 a, b, c を計算して iter の返り値のリストに (a, b, c) を追加するだけです。
それでは実行してみましょう。
*Main> primitivePythagoras 100 [(3,4,5),(5,12,13),(8,15,17),(7,24,25),(20,21,29),(9,40,41),(12,35,37)] *Main> primitivePythagoras 500 [(3,4,5),(5,12,13),(8,15,17),(7,24,25),(20,21,29),(9,40,41),(12,35,37),(11,60,61), (28,45,53),(33,56,65),(13,84,85),(16,63,65),(48,55,73),(39,80,89),(15,112,113), (36,77,85),(65,72,97),(17,144,145),(20,99,101),(60,91,109),(51,140,149),(19,180,181), (44,117,125),(88,105,137),(85,132,157),(57,176,185),(21,220,221),(24,143,145), (119,120,169),(95,168,193),(52,165,173),(104,153,185),(133,156,205),(28,195,197), (84,187,205)]
三辺の合計値が k となる全てのピタゴラス数を求めることも簡単です。つきのリストを見てください。
リスト : 合計値が k となるピタゴラス数をすべて求める
pythagoras :: Integer -> [(Integer, Integer, Integer)]
pythagoras k =
foldr check [] $ primitivePythagoras k
where check (a, b, c) xs =
let p = k `mod` (a + b + c)
q = k `div` (a + b + c)
in if p /= 0
then xs
else (q * a, q * b, q * c) : xs
primitivePythagoras k で原始ピタゴラス数を生成し、k が三辺の合計値 (a + b + c) で割り切れることを確認するだけです。とても簡単ですね。
それでは実行してみましょう。
*Main> pythagoras 12 [(3,4,5)] *Main> pythagoras 120 [(30,40,50),(20,48,52),(24,45,51)] *Main> pythagoras 240 [(60,80,100),(40,96,104),(48,90,102),(15,112,113)]
もうひとつ、クールな方法を紹介します。この方法は nineties さんのブログ ブートストラッピングでコンパイラを作る日記 Project Euler (Problem 3～10) を参考にさせていただきました。nineties さんに感謝いたします。
式 a^2 + b^2 = c^2 と a + b + c = k を使って変数 c を削除して因数分解すると、次のようになります。
a^2 + b^2 = (k - (a + b))^2 => k^2 - 2(a + b)k + (a + b)^2 - a^2 - b^2 = 0 => 2(k^2 - (a + b)k + ab) - k^2 = 0 => (k - a)(k - b) = (k^2)/2
上式より k - a と k - b は (k^2)/2 の約数であることがわかります。つまり、(k^2)/2 の約数を求め、条件 a < b < k を満たすものを探せばいいわけです。プログラムは次のようになります。
リスト : 合計値が k となるピタゴラス数をすべて求める (2)
pythagoras' :: Integer -> [(Integer, Integer, Integer)]
pythagoras' x = iter xs (reverse xs)
where xs = divisor (x * x `div` 2)
iter (y:ys) (z:zs) =
if z < y
then []
else if z >= x
then iter ys zs
else (a, b, c) : iter ys zs
where a = x - z
b = x - y
c = x - (a + b)
関数 divisor は約数をリストに昇順に格納して返します。拙作のページ Yet Another Haskell Problems 問題 54 で作成したプログラムと同じです。簡単な実行例を示します。
*Main> divisor 123456789 [1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263,123456789] *Main> divisor 100000 [1,2,4,5,8,10,16,20,25,32,40,50,80,100,125,160,200,250,400,500,625,800,1000,1250, 2000,2500,3125,4000,5000,6250,10000,12500,20000,25000,50000,100000] *Main> divisor 1000000 [1,2,4,5,8,10,16,20,25,32,40,50,64,80,100,125,160,200,250,320,400,500,625,800,1000, 1250,1600,2000,2500,3125,4000,5000,6250,8000,10000,12500,15625,20000,25000,31250, 40000,50000,62500,100000,125000,200000,250000,500000,1000000]
約数の一つを y とすると、もう一つの約数は z になります。z が y よりも大きい場合はすべての組み合わせを調べたので空リストを返します。z が x よりも大きい場合は条件を満たさないので次の組み合わせを調べます。そうでなければ条件を満たしているので、a, b, c を計算してリストに追加します。
それでは実行してみましょう。
*Main> pythagoras 1000000 [(200000,375000,425000),(218750,360000,421250)] (0.95 secs, 132363340 bytes) *Main> pythagoras' 1000000 [(200000,375000,425000),(218750,360000,421250)] (0.00 secs, 523848 bytes) 実行環境 : Windows 7, Core i7-2670QM 2.20GHz, GHCi, version 7.4.1
k が大きな値の場合、約数を高速に求めることができれば pythagoras' のほうが速くなるようです。興味のある方はいろいろ試してみてください。
実数 a の平方根 √a の値を求める場合、方程式 x2 - a = 0 を Newton (ニュートン) 法で解くことが多いと思います。方程式を f(x), その導関数を f'(x) とすると、ニュートン法は次の漸化式の値が収束するまで繰り返す方法です。
xn+1 = xn - f(xn) / f'(xn)
平方根を求める場合、導関数は f'(x) = 2x になるので、漸化式は次のようになります。
xn+1 = (xn + a / xn) / 2
参考文献 1 によると、√a より大きめの初期値から出発し、置き換え x <- (x + a / x) / 2 を減少が止まるまで繰り返すことで √a の正確な値を求めることができるそうです。
Haskell でプログラムすると、次のようになります。
リスト : 平方根を求める
sqrt' :: Double -> Double
sqrt' x
| x < 0 = error "sqrt': domain error"
| x == 0 = 0
| otherwise = iter (if x > 1 then init 1 x else 1)
where init s x
| s >= x = s
| otherwise = init (s * 2) (x / 2)
iter p = let q = (p + x / p) / 2
in if q >= p then p else iter q
局所関数 init は √x よりも大きめの初期値を求めます。たとえば、√123456 を求める場合、初期値の計算は次のようになります。
s x ------------------- 1.0 123456.0 2.0 61728.0 4.0 30864.0 8.0 15432.0 16.0 7716.0 32.0 3858.0 64.0 1929.0 128.0 964.5 256.0 482.25 512.0 241.125 √123456 = 351.363060095964
s を 2 倍、x を 1 / 2 していき、s >= x となったときの s が初期値 (512) となります。4, 16, 64, 256, ... 22n の平方根はこれだけで求めることができます。
あとは漸化式を計算して変数 q にセットし、q がひとつ前の値 p 以上になったら p を返すだけです。√123456 を求めたときの p と q の値を示します。
p q -------------------------------------- 512.0 376.5625 376.5625 352.20622925311204 352.20622925311204 351.3640693544162 351.3640693544162 351.3630600974135 351.3630600974135 351.363060095964 351.363060095964 351.363060095964 √123456 = 351.363060095964
6 回の繰り返しで √123456 を求めることができます。
整数 n の平方根の整数部分を求めることも簡単です。次のリストを見てください。
リスト : 整数 n の平方根を求める
iSqrt :: Integer -> Integer
iSqrt n
| n < 0 = error "iSqrt: domain error"
| n == 0 = 0
| otherwise = iter (init 1 n)
where init s x
| s >= x = s
| otherwise = init (s * 2) (x `div` 2)
iter p = let q = (n `div` p + p) `div` 2
in if q >= p then p else iter q
除算を / から `div` に変えただけです。√123456 と √123456789 のときの p, q の値の変化を示します。
p q (√123456 = 351) ------- 512 376 376 352 352 351 351 351 p q (√123456789 = 11111) ------------ 16384 11959 11959 11141 11141 11111 11111 11111
上記プログラム sqrt' と iSqrt では局所関数 init で初期値 s を求めていましたが、次のように引数をそのまま初期値として Newton 法を適用したほうが少し速くなりました。修正するとともにお詫び申し上げます。
リスト : 平方根を求める (修正版)
sqrt' :: Double -> Double
sqrt' x
| x < 0 = error "sqrt': domain error"
| x == 0 = 0
| otherwise = iter (if x > 1 then x else 1)
where iter p = let q = (p + x / p) / 2
in if q >= p then p else iter q
iSqrt :: Integer -> Integer
iSqrt n
| n < 0 = error "iSqrt: domain error"
| n == 0 = 0
| otherwise = iter n
where iter p = let q = (n `div` p + p) `div` 2
in if q >= p then p else iter q
もうひとつ、面白い方法を紹介しましょう。次の公式を使って平方根の整数部分を求めます。
(1) 1 + 3 + 5 + ... + (2n - 1) = n2 (2) 1 + 3 + 5 + ... + (2n - 1) = n2 < m < 1 + 3 + ... (2n - 1) + (2n + 1) = (n + 1)2
式 (1) は、奇数 1 から 2n - 1 の総和は n2 になることを表しています。式 (2) のように、整数 m の値が n2 より大きくて (n + 1)2 より小さいのであれば、m の平方根の整数部分は n であることがわかります。これは m から奇数 1, 3, 5 ... (2n - 1), (2n + 1) を順番に引き算していき、引き算できなくなった時点の (2n + 1) / 2 = n が m の平方根になります。参考文献 2 によると、この方法を「めのこ平方」と呼ぶそうです。
プログラムは次のようになります。
リスト : めのこ平方 iSqrt'' :: Integer -> Integer-> Integer iSqrt'' n m | n < m = m `div` 2 | otherwise = iSqrt'' (n - m) (m + 2)
アルゴリズムをそのままプログラムしただけなので、とくに難しいところは無いと思います。
それでは実行してみましょう。
*Main> iSqrt'' 4 1 2 *Main> iSqrt'' 16 1 4 *Main> iSqrt'' 64 1 8 *Main> iSqrt'' 80 1 8 *Main> iSqrt'' 81 1 9 *Main> iSqrt'' 82 1 9 *Main> iSqrt'' 100 1 10
この方法はとても簡単ですが、数が大きくなると時間がかかるようになります。そこで、整数を 2 桁ずつ分けて計算していくことにします。次の図を見てください。
整数 6789 を 67 と 89 に分ける 1 + 3 + ... + 15 = 82 < 67 両辺を 100 倍すると 802 < 6700 < 6789 802 = 1 + 3 + ... + 159 (= 2 * 80 - 1) 161 + 163 < (6789 - 6400 = 389) < 161 + 163 + 165
整数 6789 を 67 と 89 に分けます。最初に 67 の平方根を求めます。この場合は 8 になり、82 < 67 を満たします。次に、この式を 100 倍します。すると、802 < 6700 になり、6700 に 89 を足した 6789 も 802 より大きくなります。802 は 1 から 159 までの奇数の総和であることはすぐにわかるので、6789 - 6400 = 389 から奇数 161, 163, ... を順番に引き算していけば 6789 の平方根を求めることができます。
プログラムは次のようになります。
リスト : めのこ平方 (改良版)
iSqrt' :: Integer -> Integer
iSqrt' n
| n < 100 = iSqrt'' n 1
| otherwise = let m = 10 * iSqrt' (n `div` 100)
in iSqrt'' (n - m * m) (2 * m + 1)
iSqrt' は n の平方根の整数部分を求めます。n が 100 未満の場合は iSqrt'' で平方根を求めます。これが再帰呼び出しの停止条件になります。n が 100 以上の場合は、n の下位 2 桁を取り除いた値 (n `div` 100) の平方根を iSqrt' で求め、その値を 10 倍して変数 m にセットします。そして、iSqrt'' で n - m * m から奇数 2 * m + 1, 2 * m + 3 ... を順番に引き算していって n の平方根を求めます。
それでは実行してみましょう。
*Main> iSqrt' 6789 82 *Main> iSqrt' 123456789 11111 *Main> iSqrt' 1234567890 35136 *Main> :set +s *Main> iSqrt $ 2 ^ 100 1125899906842624 (0.00 secs, 1043332 bytes) *Main> iSqrt' $ 2 ^ 100 1125899906842624 (0.02 secs, 523936 bytes) *Main> iSqrt $ 2 ^ 10000 ・・・ 略 ・・・ (0.03 secs, 9173620 bytes) *Main> iSqrt' $ 2 ^ 10000 ・・・ 略 ・・・ (0.09 secs, 12476080 bytes) 実行環境 : Windows 7, Core i7-2670QM 2.20GHz, GHCi, version 7.4.1
iSqrt' は数が大きくなると iSqrt よりも遅くなりますが、思っていたよりも実行速度が速くて驚きました。実装がとても簡単なので、数が大きくなければ実用的に使えるかもしれません。興味のある方はいろいろ試してみてください。
ソートのお話です。拙作のページ Algorithms with Python 整列 [2] で説明した「マージソート」は、配列 buff の大きさを N とすると、大きさが N / 2 の作業用領域 work を必要としました。ここで、作業領域 work の大きさを配列 buff と同じ N にすることを考えます。この場合、最初に buff を work にコピーしておいて、再帰のたびに buff と work を交互に入れ換えることで、マージソートの実行速度を改善することができます。
なお、この方法は C++によるソート(sort)のページ 修正マージソート を参考にさせていただきました。同ページによると、『修正マージソートは、Java のクラス型のソートに採用されています。』 とのことです。有用な情報を公開されている作者様に感謝いたします。
プログラムは次のようになります。
リスト : マージソート (改良版)
def msort1(a, b, low, high):
if high - low <= 16:
insert_sort1(b, low, high)
else:
mid = (low + high) / 2
msort1(b, a, low, mid)
msort1(b, a, mid + 1, high)
i = low
j = mid + 1
k = low
while i <= mid and j <= high:
if a[i] <= a[j]:
b[k] = a[i]
i += 1
else:
b[k] = a[j]
j += 1
k += 1
while i <= mid:
b[k] = a[i]
k += 1
i += 1
while j <= high:
b[k] = a[j]
k += 1
j += 1
def merge_sort1(buff):
work = buff[:]
k = len(buff)
msort1(work, buff, 0, k - 1)
merge_sort1 は buff をコピーした配列 work を生成し、それを msort1 に渡してマージソートします。msort1 a b low high は、配列 b の区間 (low, high) を二分割してソートします。msort1 を再帰呼び出しするとき、引数 a, b を逆にすることに注意してください。二つの区間をソートしたあと、二つの区間をマージした結果は配列 a の区間 (low, high) にセットします。改良前の merge_sort では、あらかじめ buff の前半部分を work に退避していましたが、buff を work にコピーしておいて、buff と work を交互に切り替えることで、buff の前半部分を work に退避する処理が不要になります。
それでは実行結果を示します。
表 : 実行結果 (単位 : 秒)
[merge_sort]
個数 : 乱数 昇順 降順 山型
-----------------------------------
40000 : 0.189 0.107 0.167 0.140
80000 : 0.413 0.223 0.352 0.296
160000 : 0.858 0.476 0.741 0.621
320000 : 1.770 1.019 1.564 1.314
[merge_sort1 (改良版)]
個数 : 乱数 昇順 降順 山型
-----------------------------------
40000 : 0.160 0.109 0.137 0.124
80000 : 0.329 0.231 0.287 0.264
160000 : 0.729 0.511 0.624 0.577
320000 : 1.445 1.081 1.277 1.176
実行環境 : Windows 7, Core i7-2670QM 2.20GHz, Python 2.7.3
結果を見ればお分かりのように、昇順データ以外では改良版のほうが速くなりました。改良の効果は十分に出ていると思います。メモリを多く使用することになりますが、このような簡単な方法でマージソートを改良できるとは驚きました。
なお、これらの結果は M.Hiroi のコーディング、実行したマシン、プログラミング言語などの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。
# coding: shift_jis
#
# mergesort.py : マージソート
#
# Copyright (C) 2013 Makoto Hiroi
#
import time, random
# 挿入ソート
def insert_sort1(buff, low, high):
for i in xrange(low + 1, high + 1):
temp = buff[i]
j = i - 1
while j >= low and temp < buff[j]:
if temp >= buff[j]: break
buff[j + 1] = buff[j]
j -= 1
buff[j + 1] = temp
# マージソート
def msort(buff, work, low, high):
if high - low <= 16:
insert_sort1(buff, low, high)
else:
mid = (low + high) / 2
msort(buff, work, low, mid)
msort(buff, work, mid + 1, high)
#
p = 0
i = low
while i <= mid:
work[p] = buff[i]
p += 1
i += 1
i = mid + 1
j = 0
k = low
while i <= high and j < p:
if work[j] <= buff[i]:
buff[k] = work[j]
k += 1
j += 1
else:
buff[k] = buff[i]
k += 1
i += 1
while j < p:
buff[k] = work[j]
k += 1
j += 1
def merge_sort(buff):
k = len(buff)
work = [0] * (k / 2 + 1)
msort(buff, work, 0, k - 1)
# 改良版
def msort1(a, b, low, high):
if high - low <= 16:
insert_sort1(b, low, high)
else:
mid = (low + high) / 2
msort1(b, a, low, mid)
msort1(b, a, mid + 1, high)
i = low
j = mid + 1
k = low
while i <= mid and j <= high:
if a[i] <= a[j]:
b[k] = a[i]
i += 1
else:
b[k] = a[j]
j += 1
k += 1
while i <= mid:
b[k] = a[i]
k += 1
i += 1
while j <= high:
b[k] = a[j]
k += 1
j += 1
def merge_sort1(buff):
work = buff[:]
k = len(buff)
msort1(work, buff, 0, k - 1)
def test(sort, x):
# x は生成するデータの個数
# 乱数
a = [random.randint(0, 1000000) for y in xrange(x)]
# 昇順
b = range(x)
# 降順
c = range(x, 0, -1)
# 山型
d = range(x/2) + range(x/2, 0, -1)
s1 = time.clock()
sort(a)
e1 = time.clock()
print (e1 - s1)
s1 = time.clock()
sort(b)
e1 = time.clock()
print (e1 - s1)
s1 = time.clock()
sort(c)
e1 = time.clock()
print (e1 - s1)
s1 = time.clock()
sort(d)
e1 = time.clock()
print (e1 - s1)
# テスト
if __name__ == '__main__':
for n in [40000, 80000, 160000, 320000]:
test(merge_sort1, n)
クイックソートのお話です。クイックソートは、枢軸の選び方で効率が大きく左右されます。区間の中央値 [*1] を枢軸に選ぶと、区間をほぼ半分に分割することができます。この場合がいちばん効率が良く、データ数を N とすると N * log N に比例する時間でソートすることができます。
逆に、区間での最大値または最小値を枢軸に選ぶと、その要素と残りの要素の 2 つに分割にされることになります。これが最悪の場合で、分割のたびに最大値もしくは最小値を選ぶと、実行時間は要素数の 2 乗に比例することになります。つまり、単純挿入ソートと同じくらい遅いソートになります。それだけでなく、要素数が多くなるとスタックがオーバーフローする危険性もあります。
この問題は枢軸の選び方を工夫することで、完全ではありませんが回避することができます。区間の中からいくつかの要素を選び、その中の中央値を枢軸とします。たくさんの要素を選ぶとそれだけ最悪の枢軸を選ぶ危険性は少なくなりますが、値を選ぶのに時間がかかってしまいます。実際には 3 つから 5 つの要素を選んで、その中の中央値を枢軸とする場合が多いようです。
最近、M.Hiroi は Haskell の勉強で配列をソートするプログラムを作っているのですが、ネットで枢軸の選択方法を調べていたところ、9 つの要素から枢軸を選ぶ方法があることを知りました。この方法を median-of-9 と呼ぶようです。実際に 9 つの要素の中央値を求めているわけではありませんが、3 つの要素の中央値を求める方法 (median-of-3) よりも実行速度が速くなる場合があるようです。そこで、Python でプログラムを作って確かめてみました。
枢軸を選択するプログラムは次のようになります。
リスト : 枢軸の選択
# 中央値を返す
def median3(a, b, c):
if a > b:
if b > c:
return b
elif a < c:
return a
else:
return c
else:
if b < c:
return b
elif a < c:
return c
else:
return a
# 9 つの中から中央値を選ぶ
def median9(buff, low, high):
m2 = (high - low) / 2
m4 = m2 / 2
m8 = m4 / 2
a = buff[low]
b = buff[low + m8]
c = buff[low + m4]
d = buff[low + m2 - m8]
e = buff[low + m2]
f = buff[low + m2 + m8]
g = buff[high - m4]
h = buff[high - m8]
i = buff[high]
return median3(median3(a,b,c), median3(d,e,f), median3(g,h,i))
関数 median3 は引数 a, b, c の中で中央値となるものを返します。関数 median9 は区間 (low, high) から 9 つの要素を選びます。区間を (0, 1) とすると、0, 1/8, 1/4, 3/8, 1/2, 5/8, 3/4, 7/8, 1 の位置にある要素を選びます。次に、9 つの要素を 3 つのグループ (0, 1/8, 1/4), (3/18, 1/2, 5/8), (3/4, 7/8, 1) に分けて、おのおののグループの中央値を median3 で求めます。さらに、その 3 つから中央値を median3 で求め、その値が枢軸となります。今回の方法は 9 つの要素の中央値を選択しているわけではありませんが、これでも十分に効果を発揮するようです。
あとのプログラムは拙作のページ Algorithms with Python 整列 [1] 「クイックソートの改良」のプログラムとほぼ同じです。詳細は プログラムリスト をお読みください。
それでは実行結果を示します。
表 : 実行結果 (単位 : 秒)
[median-of-3]
個数 : 乱数 昇順 降順 山型
-----------------------------------
40000 : 0.115 0.049 0.053 0.189
80000 : 0.233 0.094 0.105 0.378
160000 : 0.481 0.194 0.212 0.831
320000 : 0.976 0.410 0.445 1.822
[median-of-9]
個数 : 乱数 昇順 降順 山型
-----------------------------------
40000 : 0.116 0.050 0.055 0.102
80000 : 0.241 0.107 0.112 0.210
160000 : 0.485 0.213 0.231 0.438
320000 : 0.983 0.446 0.487 0.924
実行環境 : Windows 7, Core i7-2670QM 2.20GHz, Python 2.7.3
乱数、昇順、降順のデータでは median-of-9 のほうが少し遅くなりましたが、山型データでは median-of-9 のほうが約 2 倍速くなりました。median-of-9 は少ないコストで最悪のケースを回避する優れた方法だと思います。もちろん、median-of-9 でも最悪のケースが存在するはずですが、最悪のケースに遭遇する確率は median-of-3 よりも低くなると思います。興味のある方はいろいろ試してみてください。
なお、これらの結果は M.Hiroi のコーディング、実行したマシン、プログラミング言語などの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。
# coding: shift_jis
#
# quicksort.py : クイックソートの改良
#
# Copyright (C) 2013 Makoto Hiroi
#
import time, random
# 挿入ソート
def insert_sort(buff, low, high):
for i in xrange(low + 1, high + 1):
temp = buff[i]
j = i - 1
while j >= low and temp < buff[j]:
buff[j + 1] = buff[j]
j -= 1
buff[j + 1] = temp
# 中央値を返す
def median3(a, b, c):
if a > b:
if b > c:
return b
elif a < c:
return a
else:
return c
else:
if b < c:
return b
elif a < c:
return c
else:
return a
# 9 つの中から中央値を選ぶ
def median9(buff, low, high):
m2 = (high - low) / 2
m4 = m2 / 2
m8 = m4 / 2
a = buff[low]
b = buff[low + m8]
c = buff[low + m4]
d = buff[low + m2 - m8]
e = buff[low + m2]
f = buff[low + m2 + m8]
g = buff[high - m4]
h = buff[high - m8]
i = buff[high]
return median3(median3(a,b,c), median3(d,e,f), median3(g,h,i))
# クイックソート
def quick_sort(buff, low, high):
stack = []
while True:
if high - low < 16:
insert_sort(buff, low, high)
if len(stack) == 0: break
low, high = stack.pop()
else:
pivot = median9(buff, low, high)
#pivot = median3(buff[low], buff[(low + high) / 2], buff[high])
i = low
j = high
while True:
while pivot > buff[i]: i += 1
while pivot < buff[j]: j -= 1
if i >= j: break
temp = buff[i]
buff[i] = buff[j]
buff[j] = temp
i += 1
j -= 1
if i - low > high - j:
stack.append((low, i - 1))
low = j + 1
else:
stack.append((j + 1, high))
high = i - 1
def test(sort, x):
# x は生成するデータの個数
# 乱数
a = [random.randint(0, 100000) for y in xrange(x)]
# 昇順
b = range(x)
# 降順
c = range(x, 0, -1)
# 山型
d = range(x/2) + range(x/2, 0, -1)
s1 = time.clock()
sort(a, 0, x - 1)
e1 = time.clock()
print (e1 - s1)
s1 = time.clock()
sort(b, 0, x - 1)
e1 = time.clock()
print (e1 - s1)
s1 = time.clock()
sort(c, 0, x - 1)
e1 = time.clock()
print (e1 - s1)
s1 = time.clock()
sort(d, 0, x - 1)
e1 = time.clock()
print (e1 - s1)
# テスト
if __name__ == '__main__':
for n in [40000, 80000, 160000, 320000]:
test(quick_sort, n)
1 月 5 日に出題したパズルの解答です。使用するプログラミング言語は Haskell です。今回のパズルは、演算子が + と - しかないので、数字の間に演算子を挿入して式を計算する処理は簡単にプログラムできます。ちょっと面倒なのが数字を連結する処理です。そこで、数字を連結する処理、数字の間に演算子を挿入する処理、式を計算する処理に分けてプログラムを作っていくことにします。
最初に数字を連結する処理を作ります。次のリストを見てください。
リスト : 数字の連結 concat_number :: [Int] -> [[Int]] concat_number [] = [[]] concat_number [x] = [[x]] concat_number (x:y:zs) = map (x:) (concat_number (y:zs)) ++ concat_number ((x * 10 + y):zs)
関数 concat_number は数字を格納したリストを受け取り、隣り合う数字を連結してできるパターンをすべて求めてリストに格納して返します。引数が空リストの場合は [[ ]] を返します。引数が [x] の場合は [[x]] を返します。これが再帰呼び出しの停止条件になります。
要素が 2 つ以上ある場合はリストを x : y : zs に分解して、x と y を連結しないパターンと、x と y を連結するパターンに分けて処理します。x と y を連結しない場合は x をそのまま使うことになります。リスト y : zs に concat_number を適用して、数字を連結したリストを求め、その先頭に x を追加します。この処理は map を使うと簡単ですね。x と y を連結する場合は、x * 10 + y を zs の先頭に追加し、そのリストに concat_number を適用します。あとは 2 つのリストを演算子 ++ で連結するだけです。
簡単な実行例を示します。
*Main> concat_number [1,2] [[1,2],[12]] *Main> concat_number [1,2,3] [[1,2,3],[1,23],[12,3],[123]] *Main> concat_number [1,2,3,4] [[1,2,3,4],[1,2,34],[1,23,4],[1,234],[12,3,4],[12,34],[123,4],[1234]] *Main> concat_number [1,2,3,4,5] [[1,2,3,4,5],[1,2,3,45],[1,2,34,5],[1,2,345],[1,23,4,5],[1,23,45],[1,234,5], [1,2345],[12,3,4,5],[12,3,45],[12,34,5],[12,345],[123,4,5],[123,45],[1234,5], [12345]]
次は演算子 +, - を挿入して式を生成する処理を作ります。次のリストを見てください。
リスト : 式の生成
-- 式の定義
data Expr = Val Int | Add | Sub deriving (Eq, Show)
-- 式の生成
make_expr :: [Int] -> [[Expr]]
make_expr [x] = [[Val x]]
make_expr (x:xs) = map (\zs -> (Val x):Add:zs) ys1
++ map (\zs -> (Val x):Sub:zs) ys1
where ys1 = make_expr xs
Expr で数値と演算子を定義します。Val が数値、Add が + で Sub が - です。関数 make_expr は数字を格納したリストを受け取り、数字の間に演算子を挿入するパターンをすべて求めてリストに格納して返します。
プログラムは簡単です。引数が [x] であれば、[[Val x]] を返します。そうでなければ、引数を x : xs で分解して、xs に make_expr を適用して数式を生成します。そして、その数式の先頭に map で (Val x):Add と (Val x):Sub を追加します。この処理は map を使うと簡単ですね。あとは 2 つのリストを連結するだけです。
それでは簡単な実行例を示します。
*Main> concatMap make_expr $ concat_number [1,2] [[Val 1,Add,Val 2],[Val 1,Sub,Val 2],[Val 12]] *Main> concatMap make_expr $ concat_number [1,2,3] [[Val 1,Add,Val 2,Add,Val 3], [Val 1,Add,Val 2,Sub,Val 3], [Val 1,Sub,Val 2,Add,Val 3], [Val 1,Sub,Val 2,Sub,Val 3], [Val 1,Add,Val 23], [Val 1,Sub,Val 23], [Val 12,Add,Val 3], [Val 12,Sub,Val 3], [Val 123]]
次は式を計算する処理を作ります。
リスト : 式の計算 calc_expr :: [Expr] -> Int calc_expr ((Val x):xs) = iter xs x where iter [] a = a iter (Add:(Val x):xs) a = iter xs (a + x) iter (Sub:(Val x):xs) a = iter xs (a - x)
関数 calc_expr はリストの先頭 (左側) から順番に計算していくだけです。とくに難しいところはないと思います。
簡単な実行例を示します。
*Main> map calc_expr $ concatMap make_expr $ concat_number [1,2] [3,-1,12] *Main> map calc_expr $ concatMap make_expr $ concat_number [1,2,3] [6,0,2,-4,24,-22,15,9,123] *Main> map calc_expr $ concatMap make_expr $ concat_number [1,2,3,4] [10,2,4,-4,6,-2,0,-8,37,-31,33,-35,28,20,-18,-26,235,-233,19,11,13,5,46,-22,127,119,1234]
あとは filter で指定した値になる式を取り出すだけです。プログラムは次のようになります。
リスト : 小町算の解法 komachi :: [Int] -> Int -> [[Expr]] komachi xs n = filter (\expr -> calc_expr expr == n) $ concatMap make_expr $ concat_number xs
実行例を示します。
*Main> komachi [1..9] 100 [[Val 1,Add,Val 2,Add,Val 3,Sub,Val 4,Add,Val 5,Add,Val 6,Add,Val 78,Add,Val 9], [Val 1,Add,Val 2,Add,Val 34,Sub,Val 5,Add,Val 67,Sub,Val 8,Add,Val 9], [Val 1,Add,Val 23,Sub,Val 4,Add,Val 5,Add,Val 6,Add,Val 78,Sub,Val 9], [Val 1,Add,Val 23,Sub,Val 4,Add,Val 56,Add,Val 7,Add,Val 8,Add,Val 9], [Val 12,Add,Val 3,Add,Val 4,Add,Val 5,Sub,Val 6,Sub,Val 7,Add,Val 89], [Val 12,Sub,Val 3,Sub,Val 4,Add,Val 5,Sub,Val 6,Add,Val 7,Add,Val 89], [Val 12,Add,Val 3,Sub,Val 4,Add,Val 5,Add,Val 67,Add,Val 8,Add,Val 9], [Val 123,Sub,Val 4,Sub,Val 5,Sub,Val 6,Sub,Val 7,Add,Val 8,Sub,Val 9], [Val 123,Add,Val 4,Sub,Val 5,Add,Val 67,Sub,Val 89], [Val 123,Add,Val 45,Sub,Val 67,Add,Val 8,Sub,Val 9], [Val 123,Sub,Val 45,Sub,Val 67,Add,Val 89]]
これではよくわからないので、式を文字列に変換する関数 toStr を作ります。
リスt : 式を文字列に変換
toStr :: Int -> [Expr] -> [Char]
toStr n [] = "=" ++ show n
toStr n (x:xs) =
case x of
Add -> "+"
Sub -> "-"
Val x -> show x
++ toStr n xs
プログラムは簡単なので説明は割愛します。実行例を示します。
*Main> map (toStr 100) $ komachi [1..9] 100 ["1+2+3-4+5+6+78+9=100", "1+2+34-5+67-8+9=100", "1+23-4+5+6+78-9=100", "1+23-4+56+7+8+9=100", "12+3+4+5-6-7+89=100", "12-3-4+5-6+7+89=100", "12+3-4+5+67+8+9=100", "123-4-5-6-7+8-9=100", "123+4-5+67-89=100", "123+45-67+8-9=100", "123-45-67+89=100"]
ここまでプログラムを作ると、問題を解くのは簡単です。最初の問題は次のようになります。
*Main> let a = map (komachi [1..9]) [100..999] *Main> let b = map length a *Main> maximum b 15 *Main> :m +Data.List *Main Data.List> map (+100) $ elemIndices 15 b [108,117,126]
3 桁の整数の中で、式の総数の最大値は 15 になり、その値は 108, 117, 126 の 3 通りになります。たとえば、108 になる式は次のようになります。
*Main Data.List> map (toStr 108) (a !! 8) ["1+2+3+4+5+6+78+9=108", "1+2-3+45-6+78-9=108", "1+2+34-5-6-7+89=108", "1+2+34+5+67+8-9=108", "1-2-34+56+78+9=108", "1+23+4+5+6+78-9=108", "1+23-4-5+6+78+9=108", "1+23+4+56+7+8+9=108", "12+3-4-5+6+7+89=108", "12-3+4+5-6+7+89=108", "12+3+4+5+67+8+9=108", "12+34+56+7+8-9=108", "123+4-5-6-7+8-9=108", "123-4+5-6+7-8-9=108", "123-45+6+7+8+9=108"]
解のない最小値は次のように求めることができます。
*Main Data.List> head $ map (+100) $ findIndices (==0) b 160 *Main Data.List> a !! 60 []
解が存在する最大値は次のようになります。
*Main Data.List> last $ map (+100) $ findIndices (>0) b 972 *Main Data.List> map (toStr 972) $ a !! 872 ["123+4+56+789=972"]
ちなみに、数字の並びを逆順 (9,8,7,6,5,4,3,2,1) にした場合も簡単に答えを求めることができます。
*Main Data.List> let c = map (komachi [9,8..1]) [100..999] *Main Data.List> let d = map length c *Main Data.List> maximum d 19 *Main Data.List> map (+100) $ elemIndices 19 d [102] *Main Data.List> map (toStr 102) (c !! 2) ["9+8+7+6+54-3+21=102", "9+8-7+65-4+32-1=102", "9-8+7+65-4+32+1=102", "9+8+76+5+4+3-2-1=102", "9+8+76+5+4-3+2+1=102", "9+8+76-5-4-3+21=102", "9-8+76+5-4+3+21=102", "98+7+6-5-4+3-2-1=102", "98+7+6-5-4-3+2+1=102", "98+7-6+5+4-3-2-1=102", "98+7-6+5-4+3-2+1=102", "98+7-6-5+4+3+2-1=102", "98-7+6+5+4-3-2+1=102", "98-7+6+5-4+3+2-1=102", "98-7+6-5+4+3+2+1=102", "98+7+6+5+4+3-21=102", "98-7-6-5+4-3+21=102", "98-7-6-5+43-21=102", "98+76-54+3-21=102"] *Main Data.List> head $ map (+100) $ findIndices (==0) d 194 *Main Data.List> (c !! 94) [] *Main Data.List> last $ map (+100) $ findIndices (>0) d 999 *Main Data.List> map (toStr 999) (c !! 899) ["9+8+7+654+321=999"]
パズルの世界では、1 から 9 までの数字を 1 個ずつすべて使った数字を「小町数」といいます。たとえば、123456789 とか 321654987 のような数字です。「小町算」というものもあり、たとえば 123 + 456 + 789 とか 321 * 654 + 987 のようなものです。
1 から 9 までの数字を順番に並べ、間に + と - を補って三桁の値 (100 - 999) になる式を作ることにします。100 になる式の一例を示します。
例：1 + 2 + 3 - 4 + 5 + 6 + 78 + 9 = 100
100 になる式は全部で 11 通りあります。それでは問題です。
M.Hiroi は勉強中の Haskell でプログラムを作ってみようと思います。簡単な問題ないので、興味のある方はお好みのプログラミング言語で挑戦してみてください。