完全・婚約・友愛

完全数、婚約数、友愛数についておさらい。

①完全数

ある数字nと、nのproper divisors（自分以外の約数）の和が等しいとき、nを完全数（perfect number）と呼ぶ。「博士が愛した数式」という小説にも出てきました。

例えば、6は完全数。6 = 1+2+3。

②婚約数

ある数pとqは、以下の２つの条件をともに満たすとき婚約数（betrothed number）とよばれる。準友愛数ともよばれる。

• pと、qの1以外のproper divisorsの和が等しい
• qと、pの1以外のproper divisorsの和が等しい

一番小さな婚約数は、(48, 75)。

48 = 3 + 5 + 15 + 25

75 = 2 + 3 + 4 + 6 + 8+ 12 + 16 + 24

③友愛数

ある数pとqは、以下の2つの条件をともに満たすとき友愛数(amicable number)とよばれる。

• pと、qのproper divisorsの和が等しい
• qと、pのproper divisorsの和が等しい

一番小さな友愛数は、220、284。

220 = 1 + 2 + 4 + 71 + 142

284 = 1 + 2+ 4 + 5 + 10 + 11 + 20 + 22 + 44 + 55 + 110

完全数、婚約数、友愛数を単純探索で列挙して遊んでみる。

約数はメモっておいて、婚約数、友愛数を場合は、逆引き用のテーブルを作っておけば、O(n log n)で探索できる。

#include <iostream>
#include <vector>
#include <map>

#define MAX_NUM 100000

using namespace std;

vector<int> divisors[MAX_NUM+1];
map<int, vector<int> > calcTable;

vector<int> getPerfectNumber() {
 vector<int>ret;

 for (int x = 1; x <= MAX_NUM; x++) {
   int y = 0;

   for (int i = 0; i < divisors[x].size(); i++)
     y += divisors[x][i];

   y -= x;
   if (x == y)
     ret.push_back(x);
 }
  
 return ret;
}

vector<pair<int, int> > doit(bool excludeOne) {
 map<int, vector<int> > table;

 for (int x = 1; x <= MAX_NUM; x++) {
   int y = 0;

   for (int i = 0; i < divisors[x].size(); i++)
     y += divisors[x][i];

   y -= (x+excludeOne);
   table[y].push_back(x);
 }

 vector<pair<int, int> >ret;
 for (int x = 1; x <= MAX_NUM; x++) {
   if (!table.count(x))
     continue;

   int xs = 0;
   for (int i = 0; i < divisors[x].size(); i++)
     xs += divisors[x][i];

   xs -= (x+excludeOne);
   for (int i = 0; i < table[x].size(); i++) {
     int y = table[x][i];
    
     if (y == xs && x < y)
    ret.push_back(make_pair(x, y));

   }
 }

 return ret;
}

vector<pair<int, int> > getbetrothedNumber() {
 return doit(true);
}

vector<pair<int, int> > getAmicableNumber() {
 return doit(false);
}

int main() {
 // store each number's divisors
 for (int i = 1; i <= MAX_NUM; i++) {
   for (int j = 1; j*j<=i; j++) {
     if (i % j == 0) {
    divisors[i].push_back(j);
    if (j*j != i)
      divisors[i].push_back(i/j);
     }
   }
 }

 vector<int>nums;
 vector<pair<int, int> >pairs;

 // Perfect Numbers
 nums = getPerfectNumber();
 cout << "Perfect Numbers:" << endl;
 for (int i = 0; i < nums.size(); i++)
   cout << nums[i] << endl;
 cout << endl;
 nums.clear();

 // Betrothed Numbers
 pairs = getbetrothedNumber();
 cout << "Betrothed Numbers:" << endl;
 for (int i = 0; i < pairs.size(); i++)
   cout << pairs[i].first << ", " << pairs[i].second << endl;
 cout << endl;
 pairs.clear();

 // Amicable Numbers
 pairs = getAmicableNumber();
 cout << "Amicable Numbers:" << endl;
 for (int i = 0; i < pairs.size(); i++)
   cout << pairs[i].first << ", " << pairs[i].second << endl;
 cout << endl;
 pairs.clear();

 return 0;
}

Haskellで正格評価

　遅延評価は、

・無駄な計算をしない

・無限長のデータ構造を定義できる

などの利点がある。

　「無駄な計算をしない。」については、たらいまわし関数が有名。

tarai x y z  | x <= y    = y
              | otherwise = tarai (tarai (x-1) y z) (tarai (y-1) z x) (tarai (z-1) x y)

あとは、以下のような処理に対しても遅延評価が有効的にはたらく。

(?) :: Integer -> Integer -> Integer
(?) _ 0 = 0
(?) 0 _ = 0
(?) x y = x * y

main = print $ foldl (?) 1 [100000,99999..0]

　だが、「無駄な計算をしない」=「計算が速くなる」が必ず成り立つかというとそうでもない。簡約されない式がメモリを圧迫してパフォーマンス低下の原因になることがある。たとえば、以下の式を考えてみる。

main = print $ foldl (+) 0 [1..1000000]

上の処理はスタックオーバーフローを起こす。（foldrの場合も同様。）

　なぜ、スタックオーバーが起こるのかfoldl、foldrの場合それぞれにおいて考える。

まずは、foldl、foldrを自分で定義してみる。

foldr f z []     = z
foldr f z (x:xs) = x `f` foldr f z xs

foldl f z []     = z
foldl f z (x:xs) = foldl f (f z x) xs

こうやって見てみると、foldlは命令型言語のループに近いと思う。（または、末尾再帰。）

これに対して、foldrは再帰処理のように見える。

　foldrがstack overflowを起こすのは、C++で深い再帰を書いたときにオーバーフローが起こるのと同じ原因。非末尾型の再帰のような形で式が簡約されるので、現ループの結果を出すために、次のループの結果が必要、次のループを計算するためには、次の次の結果が必要・・・。

というふうになる。結果保持しておかなければいけない値がスタックを圧迫してしまう。

　foldlの場合は、少し違う。遅延評価が関係している。foldrの場合と違って、現ループの結果はその時点で計算できるが、遅延評価のおかげで実際に必要になるまでに評価されない。

結果、thunkがヒープに書きこまれる。（簡約がはじまるまでは、遅延評価はヒープを食いつぶしている。。）実際に必要になったときに簡約が始まるが、外側から（一番深いところ）簡約が始まるため、foldlと同じことが起こってしまう。（詳しくはこのページを参照。）

　これに対してfoldl'を使うと、遅延評価をせずに、すぐにthunkをforceするため、オーバーフローは起こらない。

import Data.List

main = print $ foldl' (+) 0 [1..1000000]

foldl'は次のように定義することができる。

foldl' f z []     = z
foldl' f z (x:xs) = f z x `seq` foldl' f (f z x) xs

ここで、seq x y はxを先に簡約して、y（yはxを参照している）をreturnする関数である。

seqは以下のような書き方もできる。

f $! x = x `seq` f x

これに対して、遅延評価をするために以下のような演算子がある。（この演算子はは良く見るはず。）

f $ x  = f x

　基本的には遅延評価を使うというスタンスで、どうしても遅延評価がパフォーマンスを低下させている場合のみ、正格評価を利用するよう感じで書けばいいと思う。

遅延評価のメリットは、計算を速くすることというよりも、無限長のリストや停止しない処理を定義できることにあると思う。

Amazon Kindleで本を読もう

　Amazon Kindleすごいです。iPadみたいな専用のハードウェアを購買しないといけないかと思ってましたが、Kindle for PCという無料アプリがAmazonから配布されています。

これを使えば、パソコンで電子書籍が読めます。Kindle for PCの見た目はこんな感じです。

　Amazon.comでは、Kindle Storeというカテゴリ（下図参照）があるので、Kindle Editionの書籍を簡単に検索できます。和書と比べて、洋書は電子書籍化が速いようです。Apple StoreからiTnuesに音楽を取り入れるように、AmazonサイトからKindleに書籍を簡単に取りこむことができます。

　電子書籍は紙媒体と比べて以下の利点があります。

1. かさばらない。梅雨時は湿気の心配もいらない。
2. 輸送料が無料。（海外から取り寄せる場合も無料）
3. ワンクリックで読み始めることができる。
4. Kindle Storeの書籍は、最初の章は無料で読むことができる。（下図赤枠をクリック。）
5. ハイライトを付けることができる。
6. 一度購入すると、複数のデバイスにダウンロードできる。

　この中でも、5.の機能を推します。

この機能は、単に、自分が重要と思った場所や、おもしろいと思った場所にマークを付けることができるだけではなく、それを他の読者と共有できます。

小説であれば、ほかの読者がどういうところをおもしろいと感じたのかが分かります。専門書であれば、どこが大切なポイントになるのかが簡単に分かります。

　電子書籍に抵抗がある人も多いかと思います。私も紙媒体が好きで、電子書籍は読みにくそうと言うイメージがあったのですが、Kindleはそんなことないです。デザインもかっこよく、フォント、幅、輝度などが自分好みで調整できます。

　Kindle for PCと電子書籍のChapter. 1の試し読みは無料なので、取りあえず試してみる価値はあると思います。

Haskellで有名アルゴリズム実装（3）

今回は、クイックソート。

Haskellでのクイックソートの実装は、とても美しく、もっとも美しいコードの一つにあげる人も多いです。

qsort :: (Ord a) => [a] -> [a]
qsort []     = []
qsort (x:xs) = (qsort left) ++ [x] ++ (qsort right)
               where left  = filter (<x)  xs
                     right = filter (>=x) xs

　これで、クイックソートが出来ます。毎回リストの先頭をピボットに選択して、分割統治しています。

これだけだと、ちょっとさびしいのでクイックソートに最悪入力パターンを与えてみます。

main = print $ qsort [10000,9999..1]

　私の環境では、処理が完了するまで20秒ほどかかりました。クイックソートの理想的な場合の計算時間は長さnのリストをn/2とn/2のリストに分割して、再帰処理をするため、O(n*log n)です。が、もともとほぼ整列しているリストや、逆順にソートされているリストの場合は、nが1とn-1に分かれてしまうため、全体の計算量はO(n*n)になってしまいます。

クイックソートでは、leftとrightが均等に分かれるようなピボット選択が重要です。良く知られたピボット選択方法としては、

• ランダムな要素を選択する
• 中央の要素を選択する
• 先頭、中央、末尾の要素の中央値を選択する

などがあります。3つ目のピボット選択をHaskellで実装してみます。

selectPivot :: (Ord a) => [a] -> a
selectPivot x = let xh = head x
                    xm = x !! (div (length x) 2)
                    xl = last x
                in median [xh, xm, xl]

median :: (Ord a) => [a] -> a
median xs = let sortedXs = qsort xs
                midIndex = div (length xs) 2
            in sortedXs !! midIndex

balancedQsort :: (Ord a) => [a] -> [a]
balancedQsort [] = []
balancedQsort xs = (balancedQsort left) ++ [x] ++ (balancedQsort right)
                      where x     = selectPivot xs
                            xs'   = delete x xs
                            left  = filter  (<x) xs'
                            right = filter (>=x) xs'

上の実装だと、先ほどのデータに対しても1秒以下でソートが出来ました。恐るべしHaskell。じゃなくてピボット選択。

Haskellのページに戻る。

Church Number with SKI combinators

I wrote on my blog that lambda calculus is turing-complete the other day.

But more surprisingly you can simulate all the procedures executable on a turing machine with only three operators: S combinator, K combinator, and I combinator.

They are defined as follows respectively:

-- SKI combinator
i = \x -> x                             -- λx.x
k = \x -> \y -> x                       -- λxy.x
s = \x -> \y -> \z -> x z (y z)         -- λxyz.x z (y z)

There are other combinators, including the most famous one --Y combinator--.

I found an interesting table in which you can find a lot of combinators:

http://www.angelfire.com/tx4/cus/combinator/birds.html

Among them, I picked up a I* combinator, which seemed useful, and played around with it.

i' = s (s k)                            -- λxy.x y

Today's topic is how to represent church numbers in terms of SKI combinators.

I'll go with Haskell.

Here's the definition. I was able to define one, two, three, and four.

-- Church numeral
zero = k i
one = i
two = (s (s (k s) k)) i
three = (s (s (k s) k)) (s (s (k s) k) i)
four = (s (s (k s) k)) ((s (s (k s) k)) (s (s (k s) k) i))

I could continue it forever. But the definition above is somewhat awkward. I'll improve it later on. But for now, I'll go with it!

To assert that it really works correctly, I also implemented interfaces between lambda calculus world and standard computation world.

Here's it:

import Unsafe.Coerce

-- For assertion only
type Church = (Int -> Int) -> Int -> Int

unChurch :: Church -> Int
unChurch n = unsafeApply n (+1) (0)

unsafeApply n a b = unsafeApply' (unsafeCoerce n) a b
unsafeApply' :: Church -> (Int -> Int) -> Int -> Int
unsafeApply' n a b = n a b

main = do print $ unChurch zero
         print $ unChurch one
         print $ unChurch two
         print $ unChurch three
         print $ unChurch four
         print $ unsafeApply i' (^2) 10
         print $ unsafeApply i' (subtract 99) 10
         print $ unsafeApply ki (+40) 2

I had to use a forbiddened technique "Unsafe.Corece," since Haskell's clever type inference got in my way. For more details, let's consider the following I* combinator.

i* = s (s k)

This combinator works like λxy.x y, simply applying the function to the second argument.

And what will happen when reducing the next expression?

i* succ 0

= s (s k) succ 0

= (s k) 0 (succ 0)

= k (succ 0) (0 succ 0)

= succ 0

= 1

In the middle of the process, we've got k (succ 0) (0 succ 0). This seems OK, even though (0 succ 0) is clearly NG, because the part is not evaluated thanks to lazy evaluation.

But in Haskell, its wise type inferrence throws an error since it fails consistent type definitions.

So you have to keep GHC from inferring the type definition as for it, which can be achieved by a banned instruction "unsafeCoerce." This function seems useful when you define Y-combinator also, though Haskell has concise syntax for fixed point combinator.

Anyway you have to remmember that when using unsafeCoerce, you are responsible for its results -- your program might stop when running --. So a thorough check is required before using unsafeCoerce.

Thanks for reading.

Since I'm a newbie to Haskell, this post might include some mistakes.

In that case, please make a comment. All advices and corrections are welcomed.

Church Encoding in Haskell

　HaskellでChurch Encodingをしてみました。Church Encodingとは、すべての演算をラムダ算法でエンコードすることです。

　まずは、”数値”をラムダ式で表してみます。

これは、やってみて意味が分かりましたが、すごいことです。nという数値をf^(n) xで表すことにします。

1 -> f x

2 -> f (f x)

3 -> f (f (f x))

・・・・

n -> f^(n) x

という感じです。

fとxは何を選んでもいいですが、私たちが慣れ親しんだ数値と簡単に変換できるように

f = succ

x = 0

とします。そして、以下の2つの関数を定義します。

-- integer -> church
church n = \f -> \x -> iterate f x !! n

-- church -> integer
unChurch n = n succ 0

churchは正数をチャーチ数に変換します。unChurchはチャーチ数を対応する正数に変換します。

　これだけだと、ふーんという感じですが、このチャーチ数に対して演算を定義してあげます。例えば、チャーチ数の足し算add m nは、

　　f ^(m) x + f ^(n) x = f ^(n+m) x

となるように定義すればよいので、

add = \m -> \n -> \f -> \x -> m f (n f x)

とします。

同様に、掛け算、べき乗、インクリメントを定義します。

-- one
one = \f -> \x -> f x

-- n++
inc = \n -> \f -> \x -> f (n f x)

-- m*n
mult = \m -> \n -> \f -> \x -> m (n f) x

-- m^n
expt = \m -> \n -> n (mult m) one

　ラムダ算法の世界でも、私たちが慣れ親しんだ計算の世界と同じようにきちんと計算ができていることを確認します。

main = do let a = church 2
            b = church 10
        print $ unChurch $ add a b
        print $ unChurch $ inc a
        print $ unChurch $ mult a b
        print $ unChurch $ expt a b

　このようにラムダ算法には、”値”という概念がないにもかかわらず、チャーチ数を利用することで、基本的な計算をシュミレートすることができます。

実際、ラムダ算法はチューリング完全なので、チューリングマシン上で実行できるすべての計算処理を記述することができます。さらにいうと、すべての演算は、たった３つの演算子で記述可能です。それについてはまた今度書きます。

Haskellでメモ化する（初級）

Haskellでメモ化にチャレンジしました。

Haskellのメモ化は、C++やjavaで実装するメモ化とは違って、ひと癖もふた癖もあります。

とりあえず、基本的なメモ化ができるようになるには、2つ壁がありそう。。

1つは、参照透明性による壁。Haskellには副作用がないため、再代入ができない。つまり、メモ化する対象を更新していくというような作業ができないのです。外部変数を用意して、その外部変数に計算した値を格納していくという方法は通じないのです。。

（実は禁断のモナドを使えばできるようですが、それは最後の禁じ手だと思うので、まだ使いません。。）

2つ目は、簡約のされ方について理解が必要なことです。

ちょっと話はそれますが、このメモ化再帰を実装しているときに、改めて気がついたことがありました。

「副作用がない=外部変数を関数間で共有できない。=変数の値が勝手に変えられることはない。」

これは、大規模なプログラムを作る際に、かなり嬉しいです。というのも、参照透明性が保証されていれば、

• ある変数に何がどのタイミングで代入されるのか分からない。
• 実は関数間に順序結合があって、見た目には分かり辛い。
• 呼び出し先で勝手に値が書きかえられていた。。

などの問題から解放されるからです。これは、素晴らしい。

話を戻して、メモ化をどうやって実現するか。。更新ができなければ、毎回更新されたメモ化のコンテナを新しい変数として作成してやればいいことに気付きました。これはHaskellでは常套手段のようです。

まず、やってみたのが、これ。（例題として、フィボナッチ数列を用います。）

-- momorize status
fib n = memorized_fib n !! n

memorized_fib :: Int -> [Integer]
memorized_fib 0 = [1]
memorized_fib 1 = 1 : [1]
memorized_fib x =
   let memo = memorized_fib (x-1)
   in memo ++ [memo !! (x-1) + memo !! (x-2)]

これで、メモ化は出来ています。が、n=5000くらいになると、動きがもっさりしてきます。

Haskellのリストはランダムアクセスではなく、シーケンシャルアクセスなので、遅いのと、

簡約されない状態で値を持ち回っているのが速度を遅くしているようです。（と予想しています。）

上の例では、”値”をメモ化しましたが、Haskellでは、関数を第一級オブジェクトとして扱えるので、"関数"をメモ化することもできます。関数をメモ化するとどうなるのかやってみました。

fib n = ((memorized_fib n) !! n) n

memorized_fib :: Int -> [Int -> Integer]
memorized_fib 0 = [(\x -> 0)]
memorized_fib 1 = memorized_fib 0 ++ [(\x -> 1)]
memorized_fib n = let f  = memorized_fib (n-1)
                 v1 = (f !! (n-1)) 0
                 v2 = (f !! (n-2)) 0
             in f ++ [\x -> v1+v2]

えーっと、見てのとおり、全然だめです。一応動きますが、かなり酷いです。

引数を受け取って、その引数に対する値を返す関数をすべての引数に対してメモ化すればいいのでは？と思いましたが、きれいに書けませんでした。

そして、いろいろとネットをあさって見つけたメモ化再帰のお手本。

memorized_fib :: Int -> Integer
memorized_fib =
let fib 0 = 0
   fib 1 = 1
   fib n = memorized_fib (n-2) + memorized_fib (n-1)
in  (map fib [0 ..] !!)

無駄な処理がなく、見た目も綺麗です。速度も上の2つと比べて、高速です。

並べて比較してみると、自分のはメモ化の結果を再帰的に定義しているから余計な処理が生じているのかなと思います。あくまでもメモ化対象のものは絶対的なものとして一つ存在していて、再帰するのは実際の計算処理のみ。実際の計算と、メモ化する領域をわけて交互に連鎖させる。というかんじで作っていけばいいのでしょうか。

ただ、上のような書き方は、簡約がどのように行われるかまで把握していなければ難しい業だと思いました。Haskellでは、グラフ簡約という簡約が行われていて、一度簡約された値は、二回目以降は計算されないそうです。

ただし、すべてがそうなるわけではありません。もし、すべてがそのように処理されるのであれば、今回やっているようにわざわざ人間の手でメモ化する必要はありません。

どうやら、簡約をやり直す場合と、簡約後の値を参照しにいく、２パターンがあるようです。上のようなメモ化を書くには、”どういう記述をすれば後者のような処理になるのか。”をきちんと把握する必要がありそうです。

また、メモ化には、不動点を使ったメモ化もあるようです。データ構造にもリスト以外にtreeを使ったものがあるよう。。なぜ、そういうものが必要なのかは今の私のHaskell力では、理解できませんが、いつか分かる日が来るでしょう。。そのときに（中級）と（上級）編を書きます。