M.Hiroi's Home Page

Memorandum

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

2013 年 11月

11月3日

●循環小数

Puzzle DE Programming 循環小数 に移動しました。


2013 年 10月

10月26日

●分割数

Puzzle DE Programming 分割数 に移動しました。


10月20日

●多角数

Puzzle DE Programming 多角数 に移動しました。


10月12日

●順列の番号付け

N 個の要素の順列は N! 通りあります。この順列に 0 から N! - 1 までの番号をつけることを考えます。拙作のページ Puzzle DE Programming 幅優先探索の高速化 -- 8パズルを解く -- では、配列を使った方法を紹介しましたが、リストを使ってプログラムすることもできます。

たとえば、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 * 8!
 6 [0 1 2 3 4 5 6 7] : 8*8! + 6*7!
 7 [0 1 2 3 4 5 7] : 8*8! + 6*7! + 6*6!
 2 [0 1 2 3 4 5] : 8*8! + 6*7! + 6*6! + 2*5!
 5 [0 1 3 4 5] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4!
 4 [0 1 3 4] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3!
 3 [0 1 3] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! + 2*2!
 0 [0 1] : 8*8! + 6*7! + 6*6! + 2*5! + 4*4! + 3*3! + 2*2! + 0*1!    
 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]

正常に動作していますね。


10月6日

●ピタゴラス数

Puzzle DE Programming ピタゴラス数 に移動しました。


2013 年 9月

9月29日

●平方根の求め方

実数 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

●修正 (2013/09/30)

上記プログラム 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 よりも遅くなりますが、思っていたよりも実行速度が速くて驚きました。実装がとても簡単なので、数が大きくなければ実用的に使えるかもしれません。興味のある方はいろいろ試してみてください。

●参考文献

  1. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. 仙波一郎のページ, 『平方根計算法 (PDF)』

2013 年 5月

5月25日

●マージソートの改良

ソートのお話です。拙作のページ 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)

5月18日

●median-of-9

クイックソートのお話です。クイックソートは、枢軸の選び方で効率が大きく左右されます。区間の中央値 [*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 のコーディング、実行したマシン、プログラミング言語などの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間はかなり左右されます。興味のある方は、いろいろなデータをご自分の環境で試してみてください。

-- note --------
[*1] N 個の要素を昇順に並べたとき、中央に位置する要素 (N / 2 番目の要素) を「中央値」といいます。中央値のことを「メディアン (median) 」と呼びます。

●プログラムリスト

# 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)

2013 年 1月

1月12日

●パズルでプログラミング (解答編)

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月5日

●パズルでプログラミング

パズルの世界では、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 通りあります。それでは問題です。

  1. 式の総数が最大になる値をすべて求めてください。
  2. 解のない値で最小のものを求めてください。
  3. 解のある値で最大のものを求めてください。

M.Hiroi は勉強中の Haskell でプログラムを作ってみようと思います。簡単な問題ないので、興味のある方はお好みのプログラミング言語で挑戦してみてください。


Copyright (C) 2013 Makoto Hiroi
All rights reserved.

[ Home ]