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

Memorandum

プログラミングに関する覚え書や四方山話です。
[ Home | 2016年 1月, 3月, 7月, 9月, 11月 ]

2016 年 11 月

11月5日

●F#

F# (F Sharp) のお話です。 F# は Microsoft が開発した関数型とオブジェクト指向を融合したマルチパラダイムなプログラミング言語です。C# と同じく .NET Framework の共通中間言語 (CLI) にコンパイルされて、共通言語ランタイム (CLR) 上で実行されます。

また、F# は .NET Framework 互換の環境である Mono にも移植されているので、Mac OS や Linux などでも実行することができます。Ubuntu 系の OS であれば、次のコマンドでインストールすることができます。

$ sudo apt-get install mono-complete fsharp

●プログラムの実行

F# (Mono) にはインタプリタ (fsharpi) とコンパイラ (fsharpc) が用意されています。インタプリタは関数型言語の REPL (Read-Eval-Print-Loop) と同じで、シェルで fsharpi を起動すると、メッセージとプロンプト > が表示されます。この状態で F# のプログラムを入力して簡単に実行することができます。終了する場合は #quit;; と入力してください。

mhiroi@mhiroi-VirtualBox:~/fsharp$ fsharpi

F# Interactive for F# 4.0 (Open Source Edition)
Freely distributed under the Apache 2.0 Open Source License

For help type #help;;

> 1 + 2;;
val it : int = 3
> #quit;;

- Exit...
mhiroi@mhiroi-VirtualBox:~/fsharp$ 

コンパイラはプログラムをコンパイルして実行可能ファイル (executable file) を生成します。F# の場合、プログラムは CLI にコンパイルされるので、直接実行することはできません。mono というコマンドを使います。たとえば、Hello, World を表示するプログラムをコンパイルしてみましょう。

リスト : hello.fs

printfn "Hello World"  
mhiroi@mhiroi-VirtualBox:~/fsharp$ fsharpc hello.fs
F# Compiler for F# 4.0 (Open Source Edition)
Freely distributed under the Apache 2.0 Open Source License
mhiroi@mhiroi-VirtualBox:~/fsharp$ mono hello.exe
Hello World

ファイル名は hello.fs としました。F# の場合、ソースファイルの拡張子は .fs とするのが普通のようです。hello.fs をコンパイルすると hello.exe が生成されます。あとはシェル上で mono hello.exe を入力すると、hello.exe を実行することができます。

●簡単なベンチマーク

さて、肝心な F# (Mono) の実行速度ですが、いつものように「たらいまわし関数」を使って調べてみました。

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

open System

let rec tak x y z =
  if x <= y then z
  else
    tak (tak (x - 1) y z) (tak (y - 1) z x) (tak (z - 1) x y)

let s = Environment.TickCount
printfn "%d " (tak 22 11 0)
printfn "%d msec" (Environment.TickCount - s)
mhiroi@mhiroi-VirtualBox:~/fsharp$ fsharpc -O tarai.fs 
F# Compiler for F# 4.0 (Open Source Edition)
Freely distributed under the Apache 2.0 Open Source License
mhiroi@mhiroi-VirtualBox:~/fsharp$ mono tarai.exe
11 
2117 msec
表 : tak 22 11 0 の結果
処理系
Python3 (ver 3.4.2)98.0
Python (ver 2.7.8)89.0
Ruby (ver 2.1.2p95)50.0
Gauche (ver 0.9.3.3)33.4
PyPy (ver 2.3.1)29.0
ocamlc (ver 4.0.1)22.0
SBCL (ver 1.2.3)6.04
SML/NJ (ver 110.76)3.99
GCC -O (ver 4.9.1)2.36
GHC -O (ver 7.6.3)2.28
F# 4.0 (mono)2.12
SBCL (最適化)2.07
GCC -O2 (ver 4.9.1)1.94
ocamlopt (ver 4.0.1)1.25

F# (Mono) は Python, Ruby, PyPy とは次元の異なる速さで、ネイティブコードにコンパイルするプログラミング言語に匹敵する結果になりました。C# でも同じような結果になったので、Mono の JIT コンパイラはとても優秀だと思います。興味のある方はいろいろ試してみてください。


2016 年 9 月

9月11日

●Mono

C# (C Sharp) のお話です。 C# は Microsoft が開発したマルチパラダイムなプログラミング言語です。C# のプログラムは、同じく Microsoft が開発したフレームワーク .NET Framework の共通中間言語 (CLI) にコンパイルされて、共通言語ランタイム (CLR) 上で実行されます。

C# の実行には .NET Framework が必要になるため、Windows 以外の環境では動作しないと思われる方がいるかもしれませんが、実はそうではありません。C# や CLI は Ecma や ISO によって標準化されており、JIS 規格にも採用されています。Microsoft 以外のサードパーティーが C# コンパイラや .NET Framework 互換の環境を開発することもできるのです。その中で有名なオープンソースプロジェクトが「Mono」です。

Mono は Windows, Mac OS, Linux などマルチプラットフォームで動作する .NET Framework 互換の環境です。公式サイト からダウンロードして簡単にインストールすることができます。Ubuntu 系の OS であれば、次のコマンドでもインストールすることができます。

$ sudo apt-get install monodevelop

ただし、これは最新バージョンではないようなので、本稿では Windows 版を使うことにします。

●hello, C#!!

Mono の場合、C# のコンパイルには mcs を使います。たとえば、hello, C#!! を表示するプログラムは次のようになります。

リスト ; hello, C#!! (hello.cs)

using System;

class Test
{
  static void Main()
  {
    Console.WriteLine("hello, C#!!");
  }
}

ソースファイルの拡張子は .cs になります。msc hello.cs でソースファイルをコンパイルすると、実行ファイル hello.exe が生成されます。それを実行すると hello, C#!! と表示されます。

C>mcs hello.cs

C>dir /b hello.*
hello.cs
hello.exe

C>hello
Hello C#!

ところで、Mono で C# を学習する場合、csharp というツールを使うと便利です。csharp は Mono の配布パッケージに含まれています。csharp を起動すると、メッセージとプロンプトが表示されて、対話モードで C# が起動されます。その状態で C# の式を入力して簡単に実行することができます。

C>csharp
Mono C# Shell, type "help;" for help

Enter statements below.
csharp> quit;
C>

終了する場合は quit; (または Ctrl-D) を入力してください。

●簡単なベンチマーク

さて、肝心な C# (Mono) の実行速度ですが、いつものように「たらいまわし関数」を使って調べてみました。

リスト:たらいまわし関数 (tarai.cs)

using System;

class Tak {
  static int tak(int x, int y, int z) {
    if (x <= y) return z;
    return tak(tak(x - 1, y, z),
               tak(y - 1, z, x),
               tak(z - 1, x, y));
  }

  static void Main() {
    DateTime s = DateTime.Now;
    Console.WriteLine("{0}", tak(22, 11, 0));
    DateTime e = DateTime.Now;
    Console.WriteLine("{0}", e - s);
  }
}

それでは実行結果を示します。tak(22, 11, 0) を計算しました。他のプログラミング言語との比較を下表に示します。

表 : tak(22, 11, 0) の結果
処理系
Python (ver 2.7.3)91.9
PyPy (ver 4.0.1)17.9
SBCL (ver 1.0.55)5.85
SML/NJ (ver 110.74)3.48
GCC -O (ver 4.5.3)2.37
Julia (ver 0.4.1)2.31
C# (msc ver 4.4.2.0)2.18
SBCL (最適化)2.01
Go (ver 1.2)1.98
C# (最適化)1.98
GHC -O (ver 7.4.1)1.92
GCC -O2 (ver 4.5.3)1.89
Scala (ver 2.11.1)1.79
Java (ver 1.8.0_05)1.09
ocamlopt (ver 3.12.1)1.09

C# (Mono) は Python や PyPy とは次元の異なる速さで、ネイティブコードにコンパイルするプログラミング言語に匹敵する結果になりました。この結果には M.Hiroi も大変驚きました。Mono の JIT コンパイラはとても優秀なようです。興味のある方はいろいろ試してみてください。


2016 年 7 月

7月30日

●tagbody と go

Common Lisp のお話です。block のタグはレキシカルスコープで管理されますが、同様に tagbody と go のタグ (ジャンプ先) もレキシカルスコープで管理されます。go をラムダ式に包んで他の関数に渡すこともでき、そのラムダ式を評価するとそのタグにジャンプすることができます。

簡単な例を示しましょう。使用した処理系は SBCL です。

* (defun foo (x) (tagbody
(let ((f #'(lambda () (go exit))))
  (funcall x f) (print 1))
exit
(print 2)))

FOO
* (foo #'(lambda (f) (funcall f)))

2 
NIL
* (foo #'(lambda (f) f))

1 
2 
NIL

関数 foo の引数 x は関数で、その引数に go exit を包んだラムダ式を渡します。foo に渡す関数の中で引数を評価すると、go exit が実行されるので、tagbody のタグ exit に制御が移ります。したがって、(funcall x f) のあとの (print 1) は実行されません。引数を評価しない場合、(print 1) が実行されて、そのあとに (print 2) が実行されます。

ISLisp にも tagbody と go が用意されていて、その動作は Common Lisp と同じです。OKI-ISLisp での動作例を示します。

ISLisp>(defun print (x) (format (standard-output) "~A~%" x))
PRINT
ISLisp>(defun foo (x) (tagbody
(let ((f (lambda () (go exit))))
  (funcall x f) (print 1))
exit
(print 2)))
FOO
ISLisp>(foo (lambda (f) (funcall f)))
2
NIL
ISLisp>(foo (lambda (f) f))
1
2
NIL

7月9日

●Lisp / Scheme の繰り返し

Lisp / Scheme のお話です。繰り返しといえば、Scheme では named-let などの再帰定義、Common Lisp では dotime, dolist, do や loop でしょうか。たとえば、リストの総和を求める関数 sum-list を Scheme と Common Lisp で定義すると、次のようになります。

リスト : 総和を求める

; Scheme 版
(define (sum-list xs)
  (let loop ((xs xs) (a 0))
    (if (null? xs)
        a
        (loop (cdr xs) (+ a (car xs))))))

; Common Lisp 版
(defun sum-list (xs)
  (do ((xs xs (cdr xs))
       (a 0 (+ a (car xs))))
      ((null xs) a)))
gosh> (sum-list '(1 2 3 4 5))
15
* (sum-list '(1 2 3 4 5))

15

畳み込み (fold や reduce など) を使ったほうが簡単ですが、あえて繰り返しでプログラムしています。ここで仕様を変更して、負の要素があったら -1 を返すことにしましょう。Scheme 版は末尾再帰なので、簡単にプログラムすることができます。

リスト : 総和 (2)

; Scheme 版
(define (sum-list1 xs)
  (let loop ((xs xs) (a 0))
    (cond ((null? xs) a)
          ((negative? (car xs)) -1)
          (else (loop (cdr xs) (+ a (car xs)))))))
gosh> (sum-list1 '(1 2 3 -4 5))
-1

negative? でリストの要素 (car xs) が負かチェックして、そうであれば loop を再帰呼び出しせずに -1 を返すだけです。

Common Lisp の場合、繰り返しから脱出する return を使うと簡単です。

リスト : 総和 (2)

; Common Lisp 版
(defun sum-list1 (xs)
  (do ((xs xs (cdr xs))
       (a 0 (+ a (car xs))))
      ((null xs) a)
      (if (minusp (car xs))
          (return -1))))
* (sum-list1 '(1 2 3 -4 5))

-1

Common Lisp の場合、do 系の本体は暗黙の block (タグは nil) に囲まれていて、return や return-from nil で繰り返しから脱出することができます。リストの要素 (car xs) が負ならば、(return -1) を評価して、繰り返しから脱出します。この場合、do の返り値は return の引数 (-1) になります。

次は、リストのリストを行列とみなして、行列の要素の総和を求める関数 sum-matrix を定義しましょう。sum-matrix は負の要素を見つけたら -1 を返します。

リスト : 行列の総和

; Scheme 版
(define (sum-matrix xs)
  (let loop1 ((xs xs) (a 0))
    (if (null? xs)
        a
        (let loop2 ((ys (car xs)) (b 0))
          (cond ((null? ys)
                 (loop1 (cdr xs) (+ a b)))
                ((negative? (car ys)) -1)
                (else (loop2 (cdr ys) (+ b (car ys)))))))))

; Common Lisp 版
(defun sum-matrix (xs)
  (do ((xs xs (cdr xs))
       (a 0))
      ((null xs) a)
      (do ((ys (car xs) (cdr ys)))
          ((null ys))
          (if (minusp (car ys))
              (return-from sum-matrix -1)
            (incf a (car ys))))))
gosh> (sum-matrix '((1 2 3) (4 5 6)))
21
gosh> (sum-matrix '((1 2 3) (4 -5 6)))
-1
* (sum-matrix '((1 2 3) (4 5 6)))

21
* (sum-matrix '((1 2 3) (4 5 -6)))

-1

Scheme 版は name-let で二重ループを実現しています。負の要素を見つけたら loop1 や loop2 を再帰呼び出しせずに -1 を返すだけです。Common Lisp の場合、関数の本体は暗黙の block (タグは関数名) で囲まれているので、(return-from sum-matix -1) を評価すれば二重ループを脱出して -1 を返すことができます。

Common Lisp の場合、block のタグはレキシカルスコープです。次のように、高階関数から return-from で脱出することもできます。

リスト : 総和 (3)

; Common Lisp 版
(defun sum-list11 (xs)
  (reduce #'(lambda (a x)
              (if (minusp x)
                  (return-from sum-list11 -1)
                (+ a x)))
          xs
          :initial-value 0))
* (sum-list11 '(1 2 3 4 5 6))

21
* (sum-list11 '(-1 2 3 4 5 6))

-1
* (sum-list11 '(1 2 3 4 5 -6))

-1

畳み込みを行う関数 reduce に渡すラムダ式で、要素 x が負ならば return-from で -1 を返します。タグ sum-list11 はレキシカルスコープなので、ラムダ式の中から参照することができ、return-from でそのブロックから脱出することができます。

また、return-from tag をラムダ式に包んで他の関数に渡すこともできます。この場合、渡されたラムダ式を実行すると、tag で指定した block から脱出することができるのです。次のリストを見てください。

リスト : 行列の総和 (2)

; Common Lisp 版
(defun sum-list2 (failure xs)
  (do ((xs xs (cdr xs))
       (a 0 (+ a (car xs))))
      ((null xs) a)
      (if (minusp (car xs))
          (funcall failure))))

(defun sum-matrix1 (xs)
  (do ((xs xs (cdr xs))
       (a 0 (+ a (sum-list2 #'(lambda () (return-from sum-matrix1 -1))
                            (car xs)))))
      ((null xs) a)))
* (sum-matrix1 '((1 2 3) (4 5 6) (7 8 9)))

45
* (sum-matrix1 '((1 2 3) (4 5 6) (-7 8 9)))

-1

sum-matrix1 は行列 xs から 1 行ずつ取り出して sum-list2 に渡します。このとき、(return-from sum-matrix1 -1) を包んだラムダ式もいっしょに渡します。Common Lisp はレキシカルスコープなので、ラムダ式の中からタグ sum-matrix1 を参照することができます。次に、sum-list2 でリストの要素が負の場合、渡されたラムダ式 failure を評価します。すると、制御は sum-matrix1 に戻って -1 を返すことができます。

この動作は Scheme の call/cc による大域脱出とよく似ています。Scheme で同様のプログラムを作ると次のようになります。

リスト : 行列の総和 (2)

; Scheme 版
(define (sum-list2 failure xs)
  (let loop ((xs xs) (a 0))
    (cond ((null? xs) a)
          ((negative? (car xs)) (failure -1))
          (else (loop (cdr xs) (+ a (car xs)))))))

(define (sum-matrix1 xs)
  (call/cc
   (lambda (bk)
     (let loop ((xs xs) (a 0))
       (if (null? xs)
           a
           (loop (cdr xs) (+ a (sum-list2 bk (car xs)))))))))
gosh> (sum-matrix1 '((1 2 3) (4 5 6) (7 8 9)))
45
gosh> (sum-matrix1 '((1 2 3) (4 5 6) (7 -8 9)))
-1

call/cc で継続を取り出して変数 bk にセットし、それを sum-list2 に渡します。sum-list2 では、要素が負であれば継続 failure を評価して -1 を返します。

Common Lisp の場合、catch, throw による例外処理をサポートしているので、block とラムダ式を使って大域脱出をプログラムすることはないでしょうが、高階関数などで処理を中断させたい場合、この方法を使うことができます。それにしても、こんなことができるなんて Common Lisp は凄いですね。改めて強力なプログラミング言語だと思いました。

ところで、ISLisp は Common Lisp のサブセットなので、block tag と return-from tag でブロックから脱出することができます。ただし、Common Lisp とは違って、暗黙のブロックはありません。ご参考までに、ISLisp のプログラムと実行結果を示します。

リスト : ISLisp バージョン

(defun sum-list (xs)
  (for ((xs xs (cdr xs))
        (a 0 (+ a (car xs))))
       ((null xs) a)))

(defun sum-list1 (xs)
  (block exit
    (for ((xs xs (cdr xs))
          (a 0 (+ a (car xs))))
         ((null xs) a)
         (if (< (car xs) 0)
             (return-from exit -1)))))

; 畳み込み
(defun fold-left (f a xs)
  (for ((acc a (funcall f acc (car ys)))
        (ys xs (cdr ys)))
       ((null ys) acc)))

(defun sum-list11 (xs)
  (block exit
    (fold-left (lambda (a x)
                 (if (< x 0)
                     (return-from exit -1)
                   (+ a x)))
               0
               xs)))

(defun sum-matrix (xs)
  (block exit
    (for ((xs xs (cdr xs))
          (a 0))
         ((null xs) a)
         (for ((ys (car xs) (cdr ys)))
                  ((null ys))
                  (if (< (car ys) 0)
                      (return-from exit -1)
                    (setq a (+ a (car ys))))))))

(defun sum-list2 (failure xs)
  (for ((xs xs (cdr xs))
        (a 0 (+ a (car xs))))
       ((null xs) a)
       (if (< (car xs) 0)
           (funcall failure))))

(defun sum-matrix1 (xs)
  (block exit
    (for ((xs xs (cdr xs))
          (a 0 (+ a (sum-list2 (lambda () (return-from exit -1))
                               (car xs)))))
         ((null xs) a))))
ISLisp>(sum-list '(1 2 3 4 5))
15
ISLisp>(sum-list1 '(1 2 -3 4 5))
-1
ISLisp>(fold-left #'+ 0 '(1 2 3 4 5 6))
21
ISLisp>(sum-list11 '(1 2 3 4 5 6))
21
ISLisp>(sum-list11 '(-1 2 3 4 5 6))
-1
ISLisp>(sum-list11 '(1 2 3 4 5 -6))
-1
ISLisp>(sum-matrix '((1 2 3) (4 5 6)))
21
ISLisp>(sum-matrix '((1 2 3) (4 -5 6)))
-1
ISLisp>(sum-matrix1 '((1 2 3) (4 5 6)))
21
ISLisp>(sum-matrix1 '((1 2 3) (4 -5 6)))
-1

2016 年 3 月

3月20日

●マスターマインド (改)

M.Hiroi' Home Page で取り上げたマスターマインドは、0 から 9 までの重複しない 4 つの数字からなる隠しコードを当てるゲームでした。マスターマインドを解く場合、簡単な推測アルゴリズムを使うと、平均質問回数が 5.56 回で、質問回数の最大値は 9 回になります。

今回は数字の個数を 5 個に増やして、平均質問回数とその最大値がどうなるか、julia でプログラムを作って確かめてみました。プログラムは Julia: マスターマインドの解法 を改造すると簡単に作ることができます。説明は割愛しますので、詳細は プログラムリスト をお読みください。

結果ですが、平均質問回数が 5.99 回、質問回数の最大値は 9 で、そのときのコードは 84 通りになりました。もっと難しくなるかと思っていたので、予想外の結果にちょっと驚きました。

●プログラムリスト

#
# mastermind.jl : マスターマインドの解法
#                 (0 - 9 の数字から 5 個を選ぶ場合)
#
#                 Copyright (C) 2016 Makoto Hiroi
#

# 定数
const CSIZE = 5

# 質問したコードとその結果
type Query
    bulls::Int
    cows::Int
    code::Array{Int, 1}
end

# 0 - 9 から 5 個の数字を選ぶ順列を生成
function permutations(f, xs, n = 1)
    if n > CSIZE
        f(xs[1:CSIZE])
    else
        tmp = xs[n]
        for i in n : length(xs)
            xs[n] = xs[i]
            xs[i] = tmp
            permutations(f, xs, n + 1)
            xs[i] = xs[n]
            xs[n] = tmp
        end
    end
end

# bulls を数える
function count_bulls(xs, ys)
    c = 0
    for i in 1 : CSIZE
        if xs[i] == ys[i]; c += 1; end
    end
    c
end

# 同じ数字を数える
function count_same_number(xs, ys)
    c = 0
    for x in xs
        for y in ys
            if x == y
                c += 1
                break
            end
        end
    end
    c
end

function check(answer, xs)
    global query
    for q in query
        b = count_bulls(q.code, xs)
        c = count_same_number(q.code, xs) - b
        if b != q.bulls || c != q.cows
            return
        end
    end
    b = count_bulls(answer, xs)
    c = count_same_number(answer, xs) - b
    q = Query(b, c, xs)
    push!(query, q)
    if b == CSIZE
        throw(length(query))
    end
end

function solver()
    c = 0
    m = 0
    max_code = []
    function solver_sub(answer)
        global query
        query = Query[]
        try
            permutations(xs -> check(answer, xs), collect(0:9))
        catch e
            if m < e
                m = e
                max_code = []
            end
            if m == e
                push!(max_code, answer)
            end
            c += e
        end
    end
    permutations(solver_sub, collect(0:9))
    println(c / (10 * 9 * 8 * 7 * 6))
    println(m)
    println(max_code)
    println(length(max_code))
end

solver()
5.994246031746032
9
Any[[1,8,3,9,0],[3,9,8,0,1],[5,2,9,1,7],[5,0,6,8,3],[5,7,8,1,2],[5,8,3,7,0],
[6,5,4,1,2],[6,5,4,0,2],[6,0,1,3,9],[7,2,3,4,5],[7,3,1,4,9],[7,3,8,0,6],
[7,3,8,2,5],[7,4,0,3,5],[7,4,9,2,6],[7,5,4,0,2],[7,5,4,0,1],[7,5,9,0,3],
[7,6,8,0,3],[7,8,0,2,5],[7,8,9,6,1],[7,9,1,6,3],[8,2,1,7,9],[8,2,7,6,0],
[8,3,6,4,2],[8,4,6,0,1],[8,6,4,3,0],[8,6,5,4,3],[8,6,0,1,2],[8,7,5,4,0],
[8,7,5,0,1],[8,7,6,2,3],[8,7,6,1,2],[8,7,6,0,2],[8,7,0,6,3],[8,7,0,2,5],
[8,7,0,9,1],[8,7,9,1,2],[8,7,9,0,2],[8,0,7,6,5],[8,0,7,9,4],[8,9,1,7,2],
[9,1,0,3,8],[9,1,0,4,7],[9,2,6,0,4],[9,3,7,6,5],[9,3,8,4,0],[9,4,3,7,6],
[9,4,1,8,0],[9,4,5,0,3],[9,4,6,3,0],[9,5,3,8,7],[9,5,4,2,0],[9,5,4,8,1],
[9,5,4,8,0],[9,5,6,8,7],[9,5,7,6,8],[9,5,7,0,3],[9,6,4,0,1],[9,7,3,8,1],
[9,7,3,8,0],[9,7,4,2,5],[9,7,5,1,2],[9,7,6,1,2],[9,7,1,6,3],[9,7,0,5,3],
[9,7,0,2,5],[9,8,2,4,0],[9,8,3,5,1],[9,8,3,6,0],[9,8,3,7,1],[9,8,4,1,5],
[9,8,4,0,7],[9,8,6,1,2],[9,8,6,0,2],[9,8,7,2,5],[9,8,1,5,3],[9,8,1,0,7],
[9,8,0,5,7],[9,8,0,6,7],[9,8,0,1,7],[9,0,2,8,3],[9,0,7,5,6],[9,0,8,6,5]]
84

2016 年 1 月

1月30日

●ハミングの問題 (Hamming's Problem)

今回は「ハミングの問題」を解いてみましょう。

[ハミングの問題]

7 以上の素数で割り切れない正の整数を小さい順に N 個求めよ

参考文献 : 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991 (361 ページより引用)

7 以上の素数で割り切れない正の整数は、素因子が 2, 3, 5 しかない自然数のことです。これを「ハミング数 (Hamming Numbers)」といいます。ハミング数は素因数分解したとき、2i * 3j * 5k (i, j, k >= 0) の形式になります。たとえば、100 以下のハミング数は次のようになります。

1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 
54, 60, 64, 72, 75, 80, 81, 90, 96, 100

今回は問題を少し変えて、正の整数 n 以下のハミング数をすべて求めるプログラムを作ってみましょう。一番簡単な方法は、1 から n までの整数列を生成して、そこからハミング数を取り出していくことです。これを Python (PyPy ver 4.0.1) でプログラムすると次のようになります。

リスト : ハミングの問題

import time

def check(n):
    while n % 2 == 0: n /= 2
    while n % 3 == 0: n /= 3
    while n % 5 == 0: n /= 5
    return n == 1

def hamming(n):
    return [x for x in xrange(1, n + 1) if check(x)]

for x in xrange(2, 9):
    s = time.clock()
    print 10 ** x, len(hamming(10 ** x))
    print time.clock() - s

関数 check(n) は n がハミング数かチェックします。これは 2, 3, 5 だけで割り切れるか試しているだけです。実行結果は次のようになります。

100 34
0.00136311918243
1000 86
0.00761471403654
10000 175
0.0319288043131
100000 313
0.0208745701766
1000000 507
0.092587141429
10000000 768
0.892231480216
100000000 1105
9.03453561814

実行環境 : Windows 7, Core i7-2670QM 2.20GHz, PyPy 4.0.1

プログラムはとても簡単ですが、引数 n の値が大きくなると時間がかかるようになります。n に比べてハミング数の個数は少ないようなので、式 2i * 3j * 5k (i, j, k >= 0) を使ってハミング数を生成したほうがよさそうです。引数 n に対して i, j, k の上限値は log2 n, log3 n, log5 n で求めることができます。たとえば、100000000 の場合は次のようになります。

i : 0 - 26
j : 0 - 16
k : 0 - 11

全体で 27 * 17 * 12 = 5508 個しかありません。この中から 100000000 以下の数を選べばいいわけです。プログラムは次のようになります。

リスト : ハミングの問題 (2)

import math

def hamming2(n):
    xs2 = [2 ** x for x in xrange(0, int(math.log(n, 2)) + 1)]
    xs3 = [3 ** x for x in xrange(0, int(math.log(n, 3)) + 1)]
    xs5 = [5 ** x for x in xrange(0, int(math.log(n, 5)) + 1)]
    return sorted([x * y * z for x in xs2 for y in xs3 for z in xs5 if x * y * z <= n])

for x in xrange(8, 12):
    s = time.clock()
    print 10 ** x, len(hamming2(10 ** x))
    print time.clock() - s

2, 3, 5 のべき乗の集合を生成し、その要素を内包表記で掛け合わせて、条件を満たす数値を選択していくだけです。実行結果は次のようになりました。

100000000 1105
0.014069237311
1000000000 1530
0.00531532510766
10000000000 2053
0.0165216389339
100000000000 2683
0.0105336177752

とても速くなりましたね。このほかにも 参考文献 のプログラムが高速で、Python でジェネレータにすることも簡単です。

リスト : ハミングの問題 (3)

def hamming4():
    hs = []
    j2 = j3 = j5 = 0
    m2 = m3 = m5 = 1
    while True:
        m = min(m2, m3, m5)
        hs.append(m)
        yield m
        while m2 <= m:
            m2 = 2 * hs[j2]
            j2 += 1
        while m3 <= m:
            m3 = 3 * hs[j3]
            j3 += 1
        while m5 <= m:
            m5 = 5 * hs[j5]
            j5 += 1

for x in hamming4():
    print x,
    if x >= 100: break
1 2 3 4 5 6 8 9 10 12 15 16 18 20 24 25 27 30 32 36 40 45 48 50 54 60 64 72 75
80 81 90 96 100

このプログラムは配列 hs にハミング数をそのまま保持しているため、巨大なハミング数を求めようとするとメモリ不足になると思われます。そのような場合、不要になったハミング数を hs から削除するとうまくいくかもしれません。興味のある方はプログラムを改良してみてください。


1月9日

1 月 3 日の問題の解答です。


1月3日

●擬似完全数と不思議数

数のお話です。自然数 n の約数において、n 以外の約数の集合 xs を考えます。xs の総和が n に等しい数を「完全数」といい、xs の部分和 (部分集合の要素の和) が n に等しい数を「擬似完全数」といいます。要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題を「部分和問題」といいます。詳しい説明は拙作のページ Memorandum 2011 年 12 月 Memorandum 2012 年 1 月 の部分和問題をお読みください。

たとえば、6 の約数は {1, 2, 3, 6} で、6 以外の約数 {1, 2, 3} の総和は 6 になるので完全数です。12 の約数は {1, 2, 3, 4, 6, 12} ですが、12 以外の約数の総和は 16 になるので完全数ではありません。ところが、部分集合 {1, 2, 3, 6} の和は 12 になるので、12 は擬似完全数になります。なお、完全数は擬似完全数でもあります。

自然数 n の約数の総和が 2 * n より大きい数のことを「過剰数」といいます。擬似完全数は過剰数になります。過剰数だが擬似完全数ではない数のことを「不思議数」といいます。たとえば、70 の約数は {1, 2, 5, 7, 10, 14, 35, 70} ですが、70 を除いた約数の集合の部分和は 70 にならないので不思議数になります。

それでは問題です。n <= 10000 の範囲で、擬似完全数の個数とすべての不思議数を求めてください。


●解答

約数の求め方は拙作のページ Puzzle DE Programming 約数 で説明しています。そのとき、作成したプログラムを使いましょう。また、部分和を求めるプログラムは Memorandum 2012 年 1 月 7 日 で作成したものを使うことにすると、プログラムは次のようになります。。

リスト : 擬似完全数と不思議数

# 擬似完全数
def semiperfect_number(n):
    count = 0
    for x in xrange(2, n + 1):
        xs = divisor(x)
        if sum(xs) >= 2 * x and subset_sum(x, xs[:-1]):
            count += 1
    return count

# 不思議数
def weird_number(n):
    for x in xrange(2, n + 1):
        xs = divisor(x)
        if sum(xs) >= 2 * x and not subset_sum(x, xs[:-1]):
            print x,
    print

s = time.clock()
print semiperfect_number(10000)
print time.clock() - s
s = time.clock()
weird_number(10000)
print time.clock() -s

2 つの関数に分けましたが、一つにまとめてもかまいません。どちらの関数も約数を求める関数 divisor() を呼び出して x の約数 xs を求めます。その総和が 2 * n 以上であれば、部分和問題を解く subset_sum() を呼び出して、x と等しくなる部分集合があるかチェックします。このとき、xs から末尾の要素 (x の値) を削除することをお忘れなく。

実行結果は次のようになりました。

2485
2.01496732371
70 836 4030 5830 7192 7912 9272
1.96591159837

実行環境 : Windows 7, Core i7-2670QM 2.20GHz, PyPy 4.0.1

擬似完全数の個数は 2485 個で、不思議数の個数は 7 個になりました。実行時間は約 2 秒、divisor() だけならば 0.3 秒程度なので、部分和問題を求める subset_sum() に時間がかかるようです。部分和問題は難しいですね。興味のある方はいろいろ試してみてください。


●プログラムリスト

リスト : 擬似完全数と不思議数

import time

# 素因数分解
def factorization(n):
    def factor_sub(n, m):
        c = 0
        while n % m == 0:
            c += 1
            n /= m
        return c, n
    #
    buff = []
    c, m = factor_sub(n, 2)
    if c > 0: buff.append((2, c))
    c, m = factor_sub(m, 3)
    if c > 0: buff.append((3, c))
    x = 5
    while m >= x * x:
        c, m = factor_sub(m, x)
        if c > 0: buff.append((x, c))
        if x % 6 == 5:
            x += 2
        else:
            x += 4
    if m > 1: buff.append((m, 1))
    return buff

# 約数を求める

# p^q の約数を求める
def divisor_sub(p, q):
    a = []
    for i in xrange(0, q + 1):
        a.append(p ** i)
    return a

def divisor(n):
    xs = factorization(n)
    ys = divisor_sub(xs[0][0], xs[0][1])
    for p, q in xs[1:]:
        ys = [x * y for x in divisor_sub(p, q) for y in ys]
    return sorted(ys)

# 部分和問題
def subset_sum(n, xs):
    a = set([0])
    r = sum(xs)
    for x in xs:
        b = []
        r -= x
        for y in a:
            if x + y <= n and x + y + r >= n: b.append(x + y)
        a.update(b)
    return n in a

# 擬似完全数
def semiperfect_number(n):
    count = 0
    for x in xrange(2, n + 1):
        xs = divisor(x)
        if sum(xs) >= 2 * x and subset_sum(x, xs[:-1]):
            count += 1
    return count

# 不思議数
def weird_number(n):
    for x in xrange(2, n + 1):
        xs = divisor(x)
        if sum(xs) >= 2 * x and not subset_sum(x, xs[:-1]):
            print x,
    print

s = time.clock()
print semiperfect_number(10000)
print time.clock() - s
s = time.clock()
weird_number(10000)
print time.clock() -s

1月1日

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

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


Copyright (C) 2016 Makoto Hiroi
All rights reserved.

[ Home ]