完全・婚約・友愛

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

①完全数

ある数字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)で探索できる。

view plainprint?

1. #include <iostream>  
2. #include <vector>  
3. #include <map>  
4.   
5. #define MAX_NUM 100000  
6.   
7. using namespace std;  
8.   
9. vector<int> divisors[MAX_NUM+1];  
10. map<int, vector<int> > calcTable;  
11.   
12. vector<int> getPerfectNumber() {  
13.  vector<int>ret;  
14.   
15.  for (int x = 1; x <= MAX_NUM; x++) {  
16.    int y = 0;  
17.   
18.    for (int i = 0; i < divisors[x].size(); i++)  
19.      y += divisors[x][i];  
20.   
21.    y -= x;  
22.    if (x == y)  
23.      ret.push_back(x);  
24.  }  
25.     
26.  return ret;  
27. }  
28.   
29. vector<pair<int, int> > doit(bool excludeOne) {  
30.  map<int, vector<int> > table;  
31.   
32.  for (int x = 1; x <= MAX_NUM; x++) {  
33.    int y = 0;  
34.   
35.    for (int i = 0; i < divisors[x].size(); i++)  
36.      y += divisors[x][i];  
37.   
38.    y -= (x+excludeOne);  
39.    table[y].push_back(x);  
40.  }  
41.   
42.  vector<pair<int, int> >ret;  
43.  for (int x = 1; x <= MAX_NUM; x++) {  
44.    if (!table.count(x))  
45.      continue;  
46.   
47.    int xs = 0;  
48.    for (int i = 0; i < divisors[x].size(); i++)  
49.      xs += divisors[x][i];  
50.   
51.    xs -= (x+excludeOne);  
52.    for (int i = 0; i < table[x].size(); i++) {  
53.      int y = table[x][i];  
54.       
55.      if (y == xs && x < y)  
56.     ret.push_back(make_pair(x, y));  
57.   
58.    }  
59.  }  
60.   
61.  return ret;  
62. }  
63.   
64. vector<pair<int, int> > getbetrothedNumber() {  
65.  return doit(true);  
66. }  
67.   
68. vector<pair<int, int> > getAmicableNumber() {  
69.  return doit(false);  
70. }  
71.   
72. int main() {  
73.  // store each number's divisors  
74.  for (int i = 1; i <= MAX_NUM; i++) {  
75.    for (int j = 1; j*j<=i; j++) {  
76.      if (i % j == 0) {  
77.     divisors[i].push_back(j);  
78.     if (j*j != i)  
79.       divisors[i].push_back(i/j);  
80.      }  
81.    }  
82.  }  
83.   
84.  vector<int>nums;  
85.  vector<pair<int, int> >pairs;  
86.   
87.  // Perfect Numbers  
88.  nums = getPerfectNumber();  
89.  cout << "Perfect Numbers:" << endl;  
90.  for (int i = 0; i < nums.size(); i++)  
91.    cout << nums[i] << endl;  
92.  cout << endl;  
93.  nums.clear();  
94.   
95.  // Betrothed Numbers  
96.  pairs = getbetrothedNumber();  
97.  cout << "Betrothed Numbers:" << endl;  
98.  for (int i = 0; i < pairs.size(); i++)  
99.    cout << pairs[i].first << ", " << pairs[i].second << endl;  
100.  cout << endl;  
101.  pairs.clear();  
102.   
103.  // Amicable Numbers  
104.  pairs = getAmicableNumber();  
105.  cout << "Amicable Numbers:" << endl;  
106.  for (int i = 0; i < pairs.size(); i++)  
107.    cout << pairs[i].first << ", " << pairs[i].second << endl;  
108.  cout << endl;  
109.  pairs.clear();  
110.   
111.  return 0;  
112. }  

#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で正格評価

　遅延評価は、

・無駄な計算をしない

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

などの利点がある。

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

view plainprint?

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

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

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

view plainprint?

1. (?) :: Integer -> Integer -> Integer  
2. (?) _ 0 = 0  
3. (?) 0 _ = 0  
4. (?) x y = x * y  
5.   
6. main = print $ foldl (?) 1 [100000,99999..0]  

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

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

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

view plainprint?

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

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

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

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

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

view plainprint?

1. foldr f z []     = z  
2. foldr f z (x:xs) = x `f` foldr f z xs  
3.   
4. foldl f z []     = z  
5. foldl f z (x:xs) = foldl f (f z x) xs  

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するため、オーバーフローは起こらない。

view plainprint?

1. import Data.List  
2.   
3. main = print $ foldl' (+) 0 [1..1000000]  

import Data.List

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

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

view plainprint?

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

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は以下のような書き方もできる。

view plainprint?

1. f $! x = x `seq` f x  

f $! x = x `seq` f x

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

view plainprint?

1. f $ x  = 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でのクイックソートの実装は、とても美しく、もっとも美しいコードの一つにあげる人も多いです。

view plainprint?

1. qsort :: (Ord a) => [a] -> [a]  
2. qsort []     = []  
3. qsort (x:xs) = (qsort left) ++ [x] ++ (qsort right)  
4.                where left  = filter (<x)  xs  
5.                      right = filter (>=x) xs  

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

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

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

view plainprint?

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

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で実装してみます。

view plainprint?

1. selectPivot :: (Ord a) => [a] -> a  
2. selectPivot x = let xh = head x  
3.                     xm = x !! (div (length x) 2)  
4.                     xl = last x  
5.                 in median [xh, xm, xl]  
6.   
7. median :: (Ord a) => [a] -> a  
8. median xs = let sortedXs = qsort xs  
9.                 midIndex = div (length xs) 2  
10.             in sortedXs !! midIndex  
11.   
12. balancedQsort :: (Ord a) => [a] -> [a]  
13. balancedQsort [] = []  
14. balancedQsort xs = (balancedQsort left) ++ [x] ++ (balancedQsort right)  
15.                       where x     = selectPivot xs  
16.                             xs'   = delete x xs  
17.                             left  = filter  (<x) xs'  
18.                             right = filter (>=x) xs'  

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:

view plainprint?

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

-- 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.

view plainprint?

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

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.

view plainprint?

1. -- Church numeral  
2. zero = k i  
3. one = i  
4. two = (s (s (k s) k)) i  
5. three = (s (s (k s) k)) (s (s (k s) k) i)  
6. four = (s (s (k s) k)) ((s (s (k s) k)) (s (s (k s) k) i))  

-- 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:

view plainprint?

1. import Unsafe.Coerce  
2.   
3. -- For assertion only  
4. type Church = (Int -> Int) -> Int -> Int  
5.   
6. unChurch :: Church -> Int  
7. unChurch n = unsafeApply n (+1) (0)  
8.   
9. unsafeApply n a b = unsafeApply' (unsafeCoerce n) a b  
10. unsafeApply' :: Church -> (Int -> Int) -> Int -> Int  
11. unsafeApply' n a b = n a b  
12.   
13. main = do print $ unChurch zero  
14.          print $ unChurch one  
15.          print $ unChurch two  
16.          print $ unChurch three  
17.          print $ unChurch four  
18.          print $ unsafeApply i' (^2) 10  
19.          print $ unsafeApply i' (subtract 99) 10  
20.          print $ unsafeApply ki (+40) 2  

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つの関数を定義します。

view plainprint?

1. -- integer -> church  
2. church n = \f -> \x -> iterate f x !! n  
3.   
4. -- church -> integer  
5. unChurch n = n succ 0  

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

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

view plainprint?

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

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

とします。

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

view plainprint?

1. -- one  
2. one = \f -> \x -> f x  
3.   
4. -- n++  
5. inc = \n -> \f -> \x -> f (n f x)  
6.   
7. -- m*n  
8. mult = \m -> \n -> \f -> \x -> m (n f) x  
9.   
10. -- m^n  
11. expt = \m -> \n -> n (mult m) one  

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

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

view plainprint?

1. main = do let a = church 2  
2.             b = church 10  
3.         print $ unChurch $ add a b  
4.         print $ unChurch $ inc a  
5.         print $ unChurch $ mult a b  
6.         print $ unChurch $ expt a b  

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では常套手段のようです。

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

view plainprint?

1. -- momorize status  
2. fib n = memorized_fib n !! n  
3.   
4. memorized_fib :: Int -> [Integer]  
5. memorized_fib 0 = [1]  
6. memorized_fib 1 = 1 : [1]  
7. memorized_fib x =  
8.    let memo = memorized_fib (x-1)  
9.    in memo ++ [memo !! (x-1) + memo !! (x-2)]  

-- 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では、関数を第一級オブジェクトとして扱えるので、"関数"をメモ化することもできます。関数をメモ化するとどうなるのかやってみました。

view plainprint?

1. fib n = ((memorized_fib n) !! n) n  
2.   
3. memorized_fib :: Int -> [Int -> Integer]  
4. memorized_fib 0 = [(\x -> 0)]  
5. memorized_fib 1 = memorized_fib 0 ++ [(\x -> 1)]  
6. memorized_fib n = let f  = memorized_fib (n-1)  
7.                  v1 = (f !! (n-1)) 0  
8.                  v2 = (f !! (n-2)) 0  
9.              in f ++ [\x -> v1+v2]  

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]

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

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

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

view plainprint?

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

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力では、理解できませんが、いつか分かる日が来るでしょう。。そのときに（中級）と（上級）編を書きます。

知ってると便利なSTL(8) lower_bound、upper_bound

知ってると便利なSTLシリーズ第8弾。

今回は、○○以上の要素、○○より大きい要素の位置を返すlower_boundとupper_boundについて。

より厳密に説明すると、

lower_bound(x,y,p)は、整列済みコンテナ[x, y)の中から、与えられた値p以上となる要素の位置を返します。

upper_bound(x,y,p)は、整列済みコンテナ[x,y)の中から、与えられた値pより大きくなる要素の位置を返します。

百聞は一見に如かずということで、サンプルソース。

view plainprint?

1. #define SIZE(buff) (sizeof(buff)/sizeof(buff[0]))  
2.   
3. int main() {  
4.   int x[] = {0,1,1,2,2,2,4,5,5,10};  
5.   
6.   int p1 = (lower_bound(x, x+SIZE(x), 4) - x);  
7.   int p2 = (upper_bound(x, x+SIZE(x), 4) - x);  
8.   
9.   printf("x[%d] =%d >= %d\n", p1, x[p1], 4);  
10.   printf("x[%d] =%d > %d\n", p2, x[p2], 4);  
11.   return 0;  
12. }  

#define SIZE(buff) (sizeof(buff)/sizeof(buff[0]))

int main() {
  int x[] = {0,1,1,2,2,2,4,5,5,10};

  int p1 = (lower_bound(x, x+SIZE(x), 4) - x);
  int p2 = (upper_bound(x, x+SIZE(x), 4) - x);

  printf("x[%d] =%d >= %d\n", p1, x[p1], 4);
  printf("x[%d] =%d > %d\n", p2, x[p2], 4);
  return 0;
}

実行結果：

x[6] (=4) >= 4

x[7] (=5) > 4

さらに、逆順にソートされたコンテナから、○○以下になる要素の位置、○○より小さくなる要素の位置を取り出したいという場合もあるかと思います。そういう場合は、予めコンテナを取得または定義する段階で、マイナスを付けておくといいです。

view plainprint?

1. int main() {  
2.     // The original sequence is: {10,8,8,7,5,5,5,3,1}  
3.     int x[] = {-10,-8,-8,-7,-5,-5,-5,-3,-1};  
4.   
5.     int p1 = lower_bound(x, x+SIZE(x), -5) - x;  
6.     int p2 = upper_bound(x, x+SIZE(x), -5) - x;  
7.   
8.     printf("x[%d] (=%d) <= %d\n", p1, -x[p1], 5);  
9.     printf("x[%d] (=%d) < %d\n", p2, -x[p2], 5);  
10.   
11.     return 0;  
12. }  

int main() {
    // The original sequence is: {10,8,8,7,5,5,5,3,1}
    int x[] = {-10,-8,-8,-7,-5,-5,-5,-3,-1};

    int p1 = lower_bound(x, x+SIZE(x), -5) - x;
    int p2 = upper_bound(x, x+SIZE(x), -5) - x;

    printf("x[%d] (=%d) <= %d\n", p1, -x[p1], 5);
    printf("x[%d] (=%d) < %d\n", p2, -x[p2], 5);

    return 0;
}

これだけでは、この関数のありがたみが感じられません。何が嬉しいかというと、○○以上の要素や、○○より大きい要素の位置が対数時間で分かるということです。例えば、最長増加列の問題を高速に解きたい場合は、lower_boundが有効です。

TCO Marathon Match Round2反省記

去年に引き続き、今年もround2で敗退しそうだ。

わざわざブログに書くほどの工夫をしたわけではないが、一応やったことを書いておく。

1. 3点をランダムに選ぶ。

2. 3点の外接円をつくり、その周上にある点を選ぶ。

3. 2.で選んだ点に対して正多角形ができるかDPで探索する。

4. 凸性の検査と、重心からの距離判定を行う。

1. - 4.までをループでまわす。最初は、10角形以下は捨てるようにする。その後、この閾値を、9,8,7,..,3と下げて行く。

ちょっと工夫した点としては、

・すでにfixした点をactive集合から外して、探索から除外されるようにする。

・n角形の閾値の他に、距離閾値を設け、3点間の距離が近いものから優先探索するようにする。

・外積を利用して凸性の判定を実施することで、線形オーダーで判定する。（Testerは、O(n^2)だった。）

ぐらい。

まー、こんなもん。

上位陣は、多角形の探索を、2点と辺の数を固定して行っていた模様。確かにその方法だと、自分のやり方より明らかに速い。（自分の処理の2.3.でやりたいことが、2.だけの計算時間でできる。しかも3.は一番オーバーヘッドが大きい処理だった。。）

アイディア負け感が強いが、それ以前にに大きなミスをしていることに、コンテスト終了後に気付いた。

今回のコンテストは、相対評価のスコア。自分は、絶対ポイントの総合点が最大になるようにパラメータを設定していた。

で、こんな現象が起こってしまったと考えている。（表の数値は実際の点数と全く関係ありません。）

	case1	case2	case3	case4	case5
competitor average	100	1000	20000	300	3000
my score	80	990	22000	210	3300

これだと、total scoreではいいチューニングかもしれないが、全体としては負けてる感が。。

そもそもsmallとlargeケースで、チューニング分けろ！

てかその前に、各ケースごとに自分が作ったプログラムを評価しとけよ！

と、かなりダメダメな感じでした。

コンテスト終了後、chokudaiさんが、「上位陣と自分の各テストケースにおける優位性をしらべるために何回もsubmitしていた。」と発言しているのを見て、自分の調査の甘さにしょぼぼーーんってなりました。。

しかし、今回のmarathon matchでは、得るものが多かったのも事実です。

• オーバーヘッドの見極め方
• 最適化の仕方
• gettimeofdayの落とし穴

など今までになかった知識、方針などが身に付きました。

特に、gettimeofdayは、過去のmarathon matchでもサーバー遅いなと思いながら戦っていただけに、解決出来てかなりすっきりしました。gettimeofdayを呼ぶ回数を少なくする方法で解決したのですが、他にもいい方法があるみたいなので、それについては、いつか改めてブログに書きます。

(追記)

結果でた。まさに予想通りの現象が起きていた。平均点数が高い問題に対しては、自分の解法はかなり高い点数を取れていた。しかし、それ以外の過半数の問題に対しては平均以下のようだった。。

次回から、問題・スコアリングの特性に適したパラメータチューニングを心がけよう。

Levenshtein Distance

I recently found a classic dynamic programming problem called "Edit distance," aka,"Levenshtein distance." It's a distance between a string to another string, literally.

More specifically, how many procedure are needed to change a string A to another string B.

Allowed procedure are:

1. Insert a character into a string.

2. Delete a character from a string.

3. Change a character within a string.

For instance, consider the following situation.

You are given "kitten" as A, "sitting" as B, respectively.

0. kitten

1. sitten (change the 1st character from k to s)

2. sittin (change the 5th character from e to i)

3. sitting (insert g at the end of string)

The preceded procedure change a string "kitten" into "sitting", and also this is the minimal number of procedures needed. So the Levenshtein distance between "kitten" and "sitting" is 3.

This problem can be solved by dynamic programming. let dp[i][j] be the Levenshtein distance between substring(A, i) and substring(B, j). Here, substring(X, n) is the sub string of X between 1 and n, inclusive. (notice that n is 1-based number)

Then we come down to the next recurrence formula:

dp[i][0] = i

dp[0][j] = j

dp[i][j] = min (dp[i-1][j], dp[i][j-1], f(i, j))

where f(i, j) = 1 if A[i] == B[i]

f(i, j) = 0 if A[i] != B[i]

So, let's solve a problem on PKU. As in the problem, Levenshtein distance is helpful in dealing with DNA. I'll put my solution to it in the case my explanation above doesn't make sense to you.

view plainprint?

1. int dp[1024][1024];  
2. int main() {  
3.     int n1, n2;  
4.     string s1, s2;  
5.   
6.     while (cin >> n1 >> s1 >> n2 >> s2) {  
7.         for (int i = 0; i <= n1; i++)  
8.             dp[i][0] = i;  
9.         for (int i = 0; i <= n2; i++)  
10.             dp[0][i] = i;  
11.   
12.         for (int i = 1; i <= n1; i++) {  
13.             for (int j = 1; j <= n2; j++) {  
14.                 dp[i][j] = min(dp[i][j-1]+1, dp[i-1][j]+1);  
15.                 if (s1[i-1] == s2[j-1])  
16.                     dp[i][j] = min(dp[i][j], dp[i-1][j-1]);  
17.                 else  
18.                     dp[i][j] = min(dp[i][j], dp[i-1][j-1]+1);  
19.             }  
20.         }  
21.   
22.         cout << dp[n1][n2] << endl;  
23.     }  
24.     return 0;  
25. }  

int dp[1024][1024];
int main() {
    int n1, n2;
    string s1, s2;

    while (cin >> n1 >> s1 >> n2 >> s2) {
        for (int i = 0; i <= n1; i++)
            dp[i][0] = i;
        for (int i = 0; i <= n2; i++)
            dp[0][i] = i;

        for (int i = 1; i <= n1; i++) {
            for (int j = 1; j <= n2; j++) {
                dp[i][j] = min(dp[i][j-1]+1, dp[i-1][j]+1);
                if (s1[i-1] == s2[j-1])
                    dp[i][j] = min(dp[i][j], dp[i-1][j-1]);
                else
                    dp[i][j] = min(dp[i][j], dp[i-1][j-1]+1);
            }
        }

        cout << dp[n1][n2] << endl;
    }
    return 0;
}

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

2日目。フィボナッチ数列について。

フィボナッチ数列をどこで使うのか、何なのか？については、昔書いてるので以下のpost参照のこと。

• フィボナッチあれこれ

Haskellでフィボナッチ数列を書いてみます。いろいろ書き方があるみたいです。

まず、自分がぱっと思いついた書き方は以下のとおりです。

view plainprint?

1. fib a b = a : fib b (a+b)  

fib a b = a : fib b (a+b)

遅延評価のおかげで、フィボナッチ数列を無限番目の項まで定義できます。使うときには、実際に必要な項までを取り出せばよいです。

例えば、30番目の項までを出力させたかったら、以下のようにします。

view plainprint?

1. main = print $ take 30 $ fib 1 1  

main = print $ take 30 $ fib 1 1

次に、一般的（？）な書き方を紹介します。

view plainprint?

1. fib2 = 1 : 1 : zipWith (+) fib2 (tail fib2)  

fib2 = 1 : 1 : zipWith (+) fib2 (tail fib2)

これはクールな書き方ですね。第3項目を知るには、第1項、第2項が分かればOK。第4項目を知るには、第2項、第3項が分かればOK、・・・・。となり、一般的に、第n項を知るには第n-1項、第n-2項まで分かればOKとなっています。

この書き方も、値は必要な時点で分かっていればよいという遅延評価のメリットを活かしています。

Haskell's cool!

Haskellのページに戻る。

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

「Haskellで有名アルゴリズムを実装しよう。」企画をスタートした。三日坊主にならないように注意したい。

とりあえず、今日は簡単なやつから。

最大公約数(GCD)と最小公倍数(LCM)を求めるプログラムをHaskellで書きます。

　まず、GCD。

view plainprint?

1. gcd' :: Integer -> Integer -> Integer  
2. gcd' a 0 = a  
3. gcd' a b = gcd' b (a `mod` b)  

gcd' :: Integer -> Integer -> Integer
gcd' a 0 = a
gcd' a b = gcd' b (a `mod` b)

　実は、標準ライブラリにgcdという関数があります。自作したものはgcd'としました。

Haskellには整数値を扱う型が2つあります。IntとIntegerです。Intは32bit符号付き整数、Integerは多倍長整数です。多倍長整数が標準で使えるので便利です。

　あと、注目すべきところは、パターンマッチです。引数の”値”によって、同じ関数でも異なる定義ができます。

　次に、LCM。LCMはGCDを使って計算することができます。

view plainprint?

1. lcm' :: Integer -> Integer -> Integer  
2. lcm' a b = a * b `quot` (gcd' a b)  

lcm' :: Integer -> Integer -> Integer
lcm' a b = a * b `quot` (gcd' a b)

lcmも標準ライブラリにあるので、lcm'と関数名を付けました。

今日はここまで。Call it a day.

Haskellのページに戻る。

NとかNPとか困難とか完全とか（２）

前回、部分和問題とかナップサック問題を例にとって、NP困難とかNP完全を勉強したが、ちょっと気になることがあった。

ナップサック問題はNP-hard、部分和問題はNP-completeだが、アルゴリズマーとしては、「ん？でも動的計画で多項式時間で解けるじゃん？これってクラスPじゃないの？」って思った。
てか、NP-completeの一つであるナップサック問題が多項式時間で解けたら、他のすべてのNP-completeの問題も多項式時間で解けることになって、N=NPになるのでは？?
ん？何か勘違いしてる？

おーそうかDPで解けるのは、要素が整数の場合のみか。実数だとDPじゃ解けない。
てかこの疑似多項式時間ってのは多項式時間と何が違うの？いい情報発見。

強多項式時間問題Πに対するアルゴリズムはその実行時間のオーダーによって次のような名前が与えられる。とのこと。
1. 強多項式時間アルゴリズム (Strongly polynomial time algorithm)
時間計算量がインスタンスIの入力サイズの多項式で抑えられる

2. 弱多項式時間アルゴリズム (Weakly polynomial time algorithm)
時間計算量がインスタンスIの入力サイズ及び入力に含まれる数量のビット数の多項式で抑えられる

3. 擬多項式時間アルゴリズム (Pseudo-polynomial time algorithm)
時間計算量がインスタンスIの入力サイズ及び入力に含まれる数量の多項式で抑えられる

部分和問題の場合は、DPすると、O(ny)になるが、このyがあるから”擬多項式時間”になるわけか。(nとかyとかは前回のブログの問題定義のものです。)

てか、弱多項式時間の数量のビット数って何だ？この分かりにくい表現。要するに、数量のログオーダーって理解でいいのか？

んー、入力に含まれる数量のログオーダーっていう理解でよさそう。
ユークリッドの互除法はO((log a + log b)^2)で、これは、Weakly polynomial timeらしい。

NとかNPとか困難とか完全とか

いろいろ知識がついてきたので理解できるかなと思ったので、勉強しなおした。

まず、クラスNとNPについて。

　ともに判定問題の複雑性クラスである。つまり対象となる問題は、「Yes」、「No」で答えられる問題のみ。
判定問題の例は、ある条件を満たす整数（a,b）の組はあるか？第n番目のフィボナッチ数はm以上か？など。

よく、NPが「Not P」の略と勘違いしてしまいがちだが、これは大きな間違い。

P = deterministic Polynomial time

NP = Polynomial time

の略である。

　どっちのクラスも多項式時間で解けるが、「何を使って多項式時間で解けるか」が違う。
クラスPの問題は、決定性チューリングマシンを使って多項式時間で解ける。
一方、クラスNPの問題は、非決定性チューリングマシンを使って多項式時間で解ける。

決定性チューリングマシン(TM)とは？

　計算機を数学的に扱うために、単純化されたモデル。TMは、

・テープ

・テープの情報を読むヘッド
・機械の内部状態を記憶するメモリ

から成り立つ。"決定性"の意味は、テープから読みとった情報がxで、メモリに記憶された状態がyの場合は、毎回決まったことを実行するということ。

非決定性チューリングマシン(NTM)とは？

　テープから読みとった情報、メモリの内部状態が同じでも「毎回同じことを実行しない。次の遷移が決定されない。」

決定されないなら、どういう風に動くの？

与えられた情報から、実行しうる処理のうち、最も効率的な処理を毎回選ぶことができる。

例えば、組み合わせ最適化問題などで、今の状態から進むことのできる状態がK個ある場合、NTMは、最も早く最適解に到達できる状態へ遷移する。

ということで、分岐K、深さLの木データを探索する場合、TMの場合は、k^Lの計算が必要だが、NTMの場合は、LでOKということになる。

　ここで、注目すべきは、TMで多項式時間で解ける問題は、明らかにNTMでも多項式時間で解けるので、クラスPはクラスNPの部分集合ということである。

NP困難とNP完全

　wikipediaによると、

問題 H がNP困難のクラスに属するとは、

「NP完全問題に属するいずれかの問題 L が H へ多項式時間チューリング還元可能である」

NP困難はNP完全と同等以上に難しい。

んー、意味わからん。

ので具体例を見てみる。例として、ナップサック問題と部分和問題を取り上げる。

ナップサック問題（NP hard）
{p1, p2, ...., pn}

{c1, c2, ...., cn}

C

が与えられる。c_{i1} + c_{i2} + , ..., c_{im} <= Cとなるように{i1, i2, ..., im}を選ぶ。{p_{i1}, p_{i2}, ..., p_{im}}が最大となるような{i1, i2, ..., im}の選び方を求める。

部分和問題（NP complete）
{x1, x2, x3,.... xn}の中から任意の数字を選び、その和をyにすることができるか？

　C := y, pi := xi、ci := xiとし、ナップサック問題を解く。

もし、最大値がCならば、部分和でyが構成できる。最大値がC以下ならば、部分和でyが構成出来ない。

という手順で、部分和問題がナップサック問題へ変換できそうだ。

問題の変換は、O(n)なので、TMを用いて多項式時間でできる。

つまり、部分和問題は、ナップサック問題へ多項式時間チューリング還元可能。

　なるほど、ここまでやると「NP-hard」の定義の意味が分かった。
こうやって見ると、問題Lは問題Hの特殊形であることが分かる。問題Hが一般形の問題なので、問題Hは問題Lと同等以上に難しいという点についても納得がいく。

　ここでNP-hardはクラスNPの部分集合と誤解されがちだが、これは間違い。NP-hardの問題は、必ずしもクラスNPに属さなくてもよい。

それじゃ、次にNP-completeの定義を見てみる。

NP完全問題とは、クラスNPに属する問題でかつ、クラスNPのすべての問題から多項式時間帰着可能な問題である。

　うーん、つまりNP-completeは判定問題でなければならないと。

先ほどのNP-hardの例から分かるように、NP-completeは、すべてのNPクラスの問題の一般系と考えられる。そんな超汎用的なNP-completeは、当然NPクラスの中で一番難しいはずである。

ちょっとアバウトな議論なので、数学的に証明してみます。

NP-completeより難しい、クラスNPの問題hが存在すると仮定する。

定義より、NP-completeな問題を解くことで、問題hを解くことができる。

これは問題hがNP-copmpleteより難しいという仮定に矛盾。

よって、NP-completeはクラスNPで一番難しい。□

なんか、分かってきた気がする。

おー、そうか。

最後にwikipediaのベン図を見て、いろいろ考察する。

P!=NP予想の関係により2種類のベン図が書かれています。

・PはNPの部分集合

・NP-completeはNP-hardの部分集合(NP-completeの問題は、他の任意のNP-copmlete問題から帰着可能なため)

・クラスNPに属し、かつ、NP-hardである問題は、NP-complete

・N=NPが成り立てば、NP-complete = N(=NP)となる。

これはちょっと式を書けば分かる。（厳密な数学的証明にはなってない気はするけど。。）

NP-complete >= NP

NP-complete >= P

NP = Pより

NP-complete >= (P = NP) <= NP-complete

より、P(=NP)に含まれる問題はすべてNP-complete。　□

なんか、すっきりしたー。すっきり～～。

が、気になることが出てきた。動的計画と擬多項式時間について次のブログに書きます。

クロージャの使用例

　今日は、クロージャの使用例を書きます。昔クロージャについて書きましたが、いまいち実際どういう場面で使われるのかイメージできてませんでした。

　今回は、クロージャを使うと上手くかける例を紹介します。まず、以下のソースをみてください。（※playしてstopすると、ソースが読めます。）

クロージャの説明サンプル1 - jsdo.it - share JavaScript, HTML5 and CSS

　3つの英文があって、それを日本語訳します。上記の例では、文章が3つでしたが、文章が100個とかになるとループで書きたくなります。ループで書いてみましょう。

クロージャの説明サンプル2 - jsdo.it - share JavaScript, HTML5 and CSS

　あれ、和訳されません。何故でしょう？

alertでデバッグすると分かりますが、カウンタが3になっています。つまり、resultが返ってきたときには、すでにループを抜けて変数iは3になってるのです。

　それでは、こういうふうに書くとどうでしょうか？ループカウンタとは別にカウンタを定義して、resultを受け取った段階でインクリメントします。

クロージャの説明サンプル3 - jsdo.it - share JavaScript, HTML5 and CSS

　今度は和訳されましたが、文章の順番がおかしいです。（タイミングによっては、正常な順番で表示されるようです。）

これは、非同期処理がもたらす副作用。コールした順番にレスポンスが返ってきていないのです。

　んー困った。非同期に処理をしたいけど、ある部分では同期をとりたい。API関数を呼び出した時点でのカウンタの値のままのものを、resultが返ってきたあとも参照したい。

　ん、これは、どこかで聞いたことがある。と思ってひらめきました。

”引数以外の変数を実行時の環境ではなく、自身が定義された環境（静的スコープ）において解決することを特徴とする。”

　これって、まさにクロージャじゃないか！！

ということで、書いてみる。

クロージャの説明サンプル4 - jsdo.it - share JavaScript, HTML5 and CSS

　なるほど、そういう意味だったのか。と納得。いろんなクロージャの例をみましたが、たぶんこれが一番分かりやすいです。

Haskell勉強記（7）

foldlは、foldrで書けるらしい。

"universal property of fold"とよばれる特性を使うと変換できるらしい。読んだけど、良く分からなかった。

とりあえず、今日やったことを書きます。まずはfoldに慣れるところからはじめました。。

1) foldを自分で定義する。

view plainprint?

1. foldl' f v []     = v  
2. foldl' f v (x:xs) = foldl' f (f v x) xs  
3.   
4. foldr' f v []     = v  
5. foldr' f v (x:xs) = f x (foldr' f v xs)  
6.   
7.   
8. main = do print $ foldl' (+) 0 [1..10]  
9.        print $ foldl' (++) [] ["h", "o", "g", "e"]  
10.        print $ foldr' (:) [] ['a', 'g', 'e', 'p', 'o', 'y', 'o']  
11.        print $ foldr' (\_ n -> 1+n) 0 [1..10]  

foldl' f v []     = v
foldl' f v (x:xs) = foldl' f (f v x) xs

foldr' f v []     = v
foldr' f v (x:xs) = f x (foldr' f v xs)


main = do print $ foldl' (+) 0 [1..10]
       print $ foldl' (++) [] ["h", "o", "g", "e"]
       print $ foldr' (:) [] ['a', 'g', 'e', 'p', 'o', 'y', 'o']
       print $ foldr' (\_ n -> 1+n) 0 [1..10]

2) ラムダ式とfoldを使って標準関数を書く。

view plainprint?

1. length' = foldr (\_ n -> 1+n) 0  
2. reverse' = foldl (\xs x -> x:xs) []  
3. map' f = foldr (\x xs -> f x : xs) []  
4. filter' f = foldr (\x xs -> if f x then x : xs else xs) []  
5.   
6. main = do print $ length' [1..5]  
7.        print $ reverse' "Hello"  
8.        print $ map' (*2) [1..10]  
9.        print $ filter (<0) [0, -1, 2, -3, 4, -5]  

length' = foldr (\_ n -> 1+n) 0
reverse' = foldl (\xs x -> x:xs) []
map' f = foldr (\x xs -> f x : xs) []
filter' f = foldr (\x xs -> if f x then x : xs else xs) []

main = do print $ length' [1..5]
       print $ reverse' "Hello"
       print $ map' (*2) [1..10]
       print $ filter (<0) [0, -1, 2, -3, 4, -5]

3) セクションとfoldを使って標準関数を書く。

view plainprint?

1. sum' = foldr (+) 0  
2. product' = foldr (*) 1  
3. and' = foldr (&&) True  
4. or' = foldr (||) False  
5.   
6. main = do print $ sum' [1..10]  
7.        print $ product' [1..10]  
8.        print $ and' $ map even [2,4..10]  
9.        print $ or' $ map odd [1..10]  

sum' = foldr (+) 0
product' = foldr (*) 1
and' = foldr (&&) True
or' = foldr (||) False

main = do print $ sum' [1..10]
       print $ product' [1..10]
       print $ and' $ map even [2,4..10]
       print $ or' $ map odd [1..10]

references:

1. Writing foldl using foldr

2. A tutorial on the universality and expressiveness of fold

SRM508練習

本番は参加できなかったので、復習。

Div2 250 CandyShop

ある地点にキャンディ屋さんがある。しかしどこにあるかは覚えてはいない。i回目にキャンディ屋さんに行った時は、座標(x[i], y[i])からスタートして、R[i]回の平行移動（1回の移動では、上下左右のいづれかに1座標だけ移動可能）で辿りつけた。x, y, Rの情報を元に、キャンディ屋さんの場所を予想する。

え、、これDiv2の250??そんな簡単じゃない気が。簡単なときの500くらいのレベルあるんじゃ・・

x, yの範囲が比較的小さいので、各座標に対して矛盾がないかしらべる。具体的には、(x[i], y[i])からのマンハッタン距離がR[i]以内の場所のみを残していけばいい。

view plainprint?

1. bool ck[512][512];  
2.   
3. class CandyShop {  
4. public:  
5.    int countProbablePlaces(vector<int> X, vector<int> Y, vector<int> R) {  
6.        REP(i, 512) REP(j, 512) ck[i][j] = true;  
7.   
8.        int sz = X.size();  
9.        REP(k, sz) {  
10.            REP(i, 512) REP(j, 512)  
11.                if (ABS(X[k]+256-i) + ABS(Y[k]+256-j) > R[k]) ck[i][j] = false;  
12.        }  
13.   
14.        int cnt = 0;  
15.        REP(i, 512) REP(j, 512) if (ck[i][j]) ++cnt;  
16.   
17.        return cnt;  
18.    }  
19. };  

bool ck[512][512];

class CandyShop {
public:
   int countProbablePlaces(vector<int> X, vector<int> Y, vector<int> R) {
       REP(i, 512) REP(j, 512) ck[i][j] = true;

       int sz = X.size();
       REP(k, sz) {
           REP(i, 512) REP(j, 512)
               if (ABS(X[k]+256-i) + ABS(Y[k]+256-j) > R[k]) ck[i][j] = false;
       }

       int cnt = 0;
       REP(i, 512) REP(j, 512) if (ck[i][j]) ++cnt;

       return cnt;
   }
};

Div2 500 DivideAndShift

N個の箱からなる列がある。一番先頭からのみ物体を取ることができる。もともと先頭からM番目にある箱を取り出すには以下の操作を最低何回行えばよいか?操作は1回につき、以下のいづれかを実施できる。

1. すべての箱を後ろにシフトする。最後尾の箱は一番前にくる。

2.すべての箱を前にシフトする。一番前の箱は最後尾に行く。

3. 箱を箱の長さNの素因数pで割る。割ったあとは好きな塊を選ぶことができる。このとき、新しい長さはN/p、Mのindexは、M % (N/p) (0-based)になることに注意。

えーっと小一時間くらいやったけど、できませんでした。深さに制限をつけたDFS、それに部分メモ化再帰を追加とやりましたが、TML回避できず。これって簡単なときの900くらいはあるんじゃ・・

LayCurse様のブログを見て解決。

解法は、1.2.3.の操作はすべて可換なので、組み合わせを全探索すること。シフト操作は最後に行うものとして、3.の操作の組み合わせを探索するようです。さすが。。。

シフトと割る操作が可換なのは、モジューロ演算の交換法則からもわかります。

(a+1) % m = a % m + 1 %m (mod m)

3.の操作の順番が関係ないのは、最終的に辿り着いた箱の長さだけをみれば、indexが一意に決まることから導け出せそうです。

あと、ちょっとしたテクニックですが、M-1を再帰関数に渡すことで、0-basedのインデックスを見れるので細かい処理（こういうの m%(n/prime[i]) ? m%(n/prime[i]) : m）が必要なくなります。

view plainprint?

1. class DivideAndShift {  
2.    vector<int> prime;  
3. public:  
4.    int solve(int n, int m, int st) {  
5.        int ret = min(m, n-m);  
6.   
7.        for (int i = st; i < (int)prime.size(); i++) {  
8.            if (prime[i] > n)  
9.                break;  
10.            if (n % prime[i] == 0)  
11.                ret = min(ret, 1+solve(n/prime[i], m%(n/prime[i]), i));  
12.        }  
13.   
14.        return ret;  
15.    }  
16.   
17.    int getLeast(int N, int M) {  
18.        int n = N;  
19.   
20.        for (int i = 2; i <= N; i++) {  
21.            if (n % i == 0)  
22.                prime.push_back(i);  
23.            while (n % i == 0)  
24.                n /= i;  
25.        }  
26.   
27.        return solve(N, M-1, 0);  
28.    }  
29. };  

class DivideAndShift {
   vector<int> prime;
public:
   int solve(int n, int m, int st) {
       int ret = min(m, n-m);

       for (int i = st; i < (int)prime.size(); i++) {
           if (prime[i] > n)
               break;
           if (n % prime[i] == 0)
               ret = min(ret, 1+solve(n/prime[i], m%(n/prime[i]), i));
       }

       return ret;
   }

   int getLeast(int N, int M) {
       int n = N;

       for (int i = 2; i <= N; i++) {
           if (n % i == 0)
               prime.push_back(i);
           while (n % i == 0)
               n /= i;
       }

       return solve(N, M-1, 0);
   }
};

この問題は、計算量を見積るのが難しい気がするのですが、どうなんでしょう？

素因数の数の上限が20くらいなので、2^20??　であってますかね？

Yコンビネータ実装しなイカ？

Yコンビネータについて勉強しなイカ？

ということで、勉強してみました。

　まず、Yコンビネータとは何か？それを知るには、少し導入が必要になります。

数学的な話になりますが、不動点というものがあります。不動点とは、ある写像に対して自分自身に写される点のことです。例えば、f(x) = x^2という関数に対して、x = 0、x = 1という点は不動点です。（f(0) = 0、f(1) = 1だからです。）

　さて、上の説明では一階関数の不動点を例にとりました。同じように高階関数に対しても不動点が存在します。高階関数fと、別の関数pに対してf(p) = pが成り立つとき、関数pを関数fの不動点といいます。

関数gを入力に取り、関数gの不動点を返す高階関数を不動点コンビネータと呼びます。

　Y(g) = g(Y(g))

となるようなYが不動点コンビネータです。

そして、代表的な不動点コンビネータが、名字も名前も有名なかのハスケル・カリーさんが発見したYコンビネータ

　Y = λf.(λx.f (x x)) (λx.f (x x))　

です。

いきなりλが出てきて謎な感じですが、このλはλ算法と呼ばれるものです。λ算法では、関数x+2をλx.x+2のように書きます。また、λ算法は、カリー化された表現なので、x-yをλxy.x-yのように書きます。

Yコンビネータに関数gを適用させてみます。

　Yg = (λx. g(x x)) (λx. g(x x))

　　= g ((λx. g(x x)) (λx. g(x x)))

　　 = g (Yg)

となります。1行目から2行目にかけては、β簡約が行われています。

　では、HaskellでYコンビネータを実装してみます。

しかし、Haskellでは遅延評価のおかげで、不動点コンビネータがそのまま書けてしまえます。Haskell神です。みんなもHaskell使ってみなイカ？

view plainprint?

1. fix :: (t -> t) -> t  
2. fix f = f (fix f)  

fix :: (t -> t) -> t
fix f = f (fix f)

のように定義します。この定義により、明らかにfixは関数fの不動点コンビネータになっています。

プログラミングの世界では、不動点コンビネータを利用することにより、無名再帰ができます。名前のないはずの関数が、名前のないはずの自分自身を呼び出すというちょっと奇妙なトリックですが、以下のようにできます。

view plainprint?

1. main = print $  
2. 　　　　(fix (\f x -> if x == 0 then 1 else x * f (x-1))) 10  

main = print $
　　　　(fix (\f x -> if x == 0 then 1 else x * f (x-1))) 10

上のmainアクションでは、10の階乗を出力しています。

ここでポイントとなるのは、無名関数の書き方です。

普通、階乗を書くとしたら、

view plainprint?

1. fact = \x -> if x == 0 then 1 else x * fact (x-1)  

fact = \x -> if x == 0 then 1 else x * fact (x-1)

ですが、不動点コンビータで使用するためには、

view plainprint?

1. fact = \f x -> if x == 0 then 1 else x * f (x-1)  

fact = \f x -> if x == 0 then 1 else x * f (x-1)

のように書きます。

関数がどのように簡約されていくかをみれば、なぜこのような書き方をするのかが分かります。

(fix fact) 10

⇒ fact (fix fact) 10

⇒ 10 * ((fix fact) 9)

⇒ 10 * (fact (fix fact) 9)

⇒ 10 * (9 * ((fix fact) 8))

・・・・

で10の階乗が計算できるわけです。

　しかし、この書き方だと、fix自身が名前付きの再帰となっていて、ちょっとインチキくさいです。

ここの記事にも、同じようなことが書いてありました。

「ラムダ計算がチューリング完全であることを示すのが目的であれば、この定義はインチキです。」

うわっ、なんかすごい学術的な言葉が出てきた。以下wikipediaから引用。

"チューリング完全とは、あるプログラミング言語がチーリング機械と同じ計算能力を持つときその言語はチューリング完全であるという。"

な、なるほど。ここで言いたいのは、ラムダ計算には名前付き再帰がないにも関わらず、名前付きの再帰を使っていては、チューリング完全を満たさない。ということでしょう。。先のブログでは、無名関数による実装が紹介されていました。

自分も、もう少しHaskell力がついたらチャレンジしたいです。