Page List

Search on the blog

2011年7月31日日曜日

カーマイケル数の列挙

カーマイケル数を列挙してみる。カーマイケル数といえば、フェルマーテストの偽陽性を引き起こすクセモノ。あらかじめカーマイケル数を列挙しておけば、フェルマーテストと併用することで、正確な素数判定ができる。

いろいろな関数を使うことになった。
  • 最大公約数
  • 最小公倍数
  • 素因数分解
  • 繰り返し二乗法
など基本的な数学系処理の練習にいいかも。。

さて、問題のカーマイケル数の部分だが、カーマイケル関数を使って、


となる最小のmを求める。このmをλ(n)と表記する。ここで、aとnはcoprimeである。カーマイケル数は、


となるnなので、(n-1)がλ(n)で割り切れるとき、nはカーマイケル数の候補となる。さらに、ここから素数を除く必要がある。素数の場合は、φ(n) = n-1となるので、それを除けばいい。逆に、カーマイケル数は、合成数なので、totient functionはn-1より小さくなるため、必ずφ(n) != n-1となり、この除去で偽陰性となることはない。
map<int, int>factorize(int x) {
map<int, int> ret;

for (int i = 2; i*i <= x; i++) {
while (x % i == 0) {
ret[i]++;
x /= i;
}
}

if (x != 1)
ret[x] = 1;

return ret;
}

int mpow(int x, int p) {
int ret = 1;

while (p) {
if (p & 1)
ret = ret * x;
x = x * x;
p >>= 1;
}

return ret;
}

int gcd(int a, int b) {
if (b)
return gcd(b, a % b);
return a;
}

int lcm(int a, int b) {
return (long long int) a * b / gcd(a, b);
}

int carmichaelFunction(int x) {
map<int, int> fact = factorize(x);

int ret = 1;
for (map<int, int>::iterator itr = fact.begin(); itr != fact.end(); itr++) {
int p = itr->first;
int k = itr->second;

int lambda = (p > 2 || k < 3) ? mpow(p, k-1) * (p-1) : mpow(p, k-2);
ret = lcm(ret, lambda);
}

return ret;
}

int eulerFunction(int x) {
map<int, int> fact = factorize(x);

double ret = x;
for (map<int, int>::iterator itr = fact.begin(); itr != fact.end(); itr++)
ret *= (1 - 1./itr->first);

return ret + 0.5;
}

int main() {

for (int i = 2; i < 100000; i++) {
int lambda = carmichaelFunction(i);
int phi = eulerFunction(i);

if ((i-1) % lambda == 0 && phi != i-1) {
printf("%d\n", i);
}
}

return 0;
}

カーマイケル数は、561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361...と続く。

2011年7月29日金曜日

拡張ユークリッド互除法について

 蟻本を見ていると、以下のような文章がでてくる。

「(省略) これは、a * x + b * y = 1 となる(x, y)を求める問題に他ならない。この式から明らかなようにgcd (a, b) != 1のときは、整数解(x, y)は存在しません。」

 自分にとっては明らかじゃなかったので、ちょっと考えてみた。
gcd(a, b) != 1のとき、上記の条件を満たす(x, y)が存在しないことを背理法で示す。

gcd(a, b) = d, d > 1とすると、
a * x + b * y = 0 (mod d)
である。
しかし、
1 != 0 (mod d)
であるため、矛盾。よって、gcd(a, b) != 1のときは上記の条件を満たす(x, y)は存在しない。□

あら、確かに。自明レベルの内容だった。。

 逆に、gcd(a, b) = 1 のときは、上記の条件を満たす(x, y)が必ず存在することを証明してみる。
a * x + b * y = 1より
x = (1 - b * y) / a となる(x, y)を求めればよい。
つまり、1- b * y = 0 (mod a)となる任意の整数yが存在すればよいことになる。
移行すると、
b* y = 1 (mod a)
となるが、gdc(a, b) = 1 からaとbがrelatively primeなので、bの逆元yは存在する。
よって、gcd(a, b) = 1のときは、解(x, y)が存在する。□

ふむふむ。なるほど。。

 証明が出来てもプログラムが書けないことには意味がない。ということで書いてみる。
int extGcd(int a, int b, int &x, int &y) {
int d = a;

if (b) {
d = extGcd(b, a%b, y, x);
y -= a / b * x;
} else {
x = 1;
y = 0;
}

return d;
}
蟻本では、b == 0のとき、y = 0としているが、ここは任意の整数でいい。(逆元は一つだけじゃないので。)たとえば、(a, b) = (11, 13)として、a*x + b * y = 1となる(x, y)を複数個求めてみる。

ちょっとだけ、extGcd()を書きかえて、b=0のときのyの値を選べるようにしてやる。(本当は、クロージャみたいなので書くといいのかも。Haskellでやればよかったかな・・)
int extGcd(int a, int b, int y0, int &x, int &y) {
int d = a;

if (b) {
d = extGcd(b, a%b, y0, y, x);
y -= a / b * x;
} else {
x = 1;
y = y0;
}

return d;
}

int main () {
int x, y;
int a = 11;
int b = 13;

for (int i = 0; i < 10; i++) {
int d = extGcd(a, b, i, x, y);
printf("%d * %d + %d * %d = %d\n", a, x, b, y, (a*x+b*y));
}
return 0;
}
で、(x, y)の組が10個求まりました。

2011年7月28日木曜日

Sushida Killer

I've recently developed an application named "Sushida Killer," which automatically plays "Sushida," the most famous typing game here in Japan -- For more details, play sushida here --.

This application consists of three parts:
1. Reading problems
2. Recognizing characters in movie
3. Similate key typing

The first part is reading problems.
On first sight, you will notice that a dish of sushi, which corresponds to one problem, emerges every 5 second. So it's not a bad idea to automatically take a screen shot of flash video with a span of 5 seconds. How to get a screen shot automatically? I used the next procedures:

1. Simulate key typing "Alt + PrintScrenn"
2. Access to the clip board with win32 API.
3. Get the Bitmap infomation from the clip board
4. Save it to a bmp file.

The second part is the key part of Sushida Killer.
To read only the charactor on screen, I first created some mimic cpp library that deals with image processing.
With it, I cut off a rectangle, in which a problem text comes up, then changed RGB 3-demension vector infomation of each pixel into 0/1 data using a threshold.
So far, you can get data something like this:
I had to collect many sample data like above to train the classifier, which was done semi-automatically.
After training with enough amounts of learning data sets, the next step is classifying an unknown data to a string. Fortunately I didn't have to use a learning machine based on statistical theory like an SVM and an Bayes' machine, since the characters were written by Flash app, not by a human, so they were neat and tidy enough to be understood by computer.

The key point is speed. How to classify data in less than 5 second?
There are about 1000 sample data, each has a 10,000-pixel relevant rectangle in it.
It seems tough to finish classifying the image data to the corresponding string in 5 seconds, considering that you also have to get a screen shot of the flash movie, and save it to the local file, then classify it into one of about 1000 groups.
So I mapped each learning data's information to a hash value between 0 and 10^10, inclusive.
I thought a collision would never happen because (the number of data / hash value size) was 10^-7. And I was right!! I succeed to finish the process less than 5 second. Actually, it doesn't take more than 1 second.
I preliminarily calculated each data's hash value, and dispatched them to a file. So when playing Sushida, program can load the data in memory, which doesn't take time at all.

The third part is key simulating.
It's not so difficult to simulate key typing if you know about keybd_event() of win32 API.
But getting the window handle of Sushida was a little challenging. The flash app runs in a browser, so it's impossible to get the application window directly.
I had to implement a program that gives you a window handle of the application on which the mouse cursor is located. Thanks to it, I was able to get the handle of Sushida, which was the children of children of children of children of my web browser!!! So I traversed window handle to reach the flash app.

Most of my effort was put into collecting the data. That were tedious tasks, though I did it semi-automatically.
And also selecting parameters was not so interesting, but it was a very important job.
What the size of rectangle?
How long the program have to sleep to wait for the screen shot and file writing to end?
How long the program have to sleep to wait for the key simulation to end?

Finally, I thank the great developer who created sushida for giving me this good challenge.
(Ah, he/she never gave me that kind of challenge, I just think I was given it..)
And I hereby promise that I will never use this Sushida killer to be ranked in the formal site, i.e., this program is for local use only.

Last but not least, here's a demo of Sushida Killer!


2011年7月26日火曜日

無理数を有理数で近似

無理数xを、有理数N/Dで近似するというのは、よくあるテーマ。
現実的には、適当に有理数の分子、分母に上限を付ける。例えば、1 <= N, D <= Lのような範囲のN, Dに対して、|x - N/D|が最小となるようなN, Dを求める問題を考えることになる。

何も考えないで実装すると、N, Dを全探索するはめになり、O(L^2)の擬多項式時間計算量が必要となる。アルゴリズムをやったことがある人なら、Dを全部まわして、Nを二分探索するという方法を思いつくだろう。これだと、O(L * log L)で L = 1,000,000程度でも一瞬で答えが出る。

実は二分探索よりも、もっと面白い方法がある。
N/D > x なら 探索点を N/(D+1)に移動、そうでなければ探索点を(N+1)/Dに移動。これを逐次繰り返していき、最良点を保持すれば答えが出る。

以下、サンプルソース。無理数じゃないけど、x = 100007/100009として解いてみる。
#define MAX_VAL 1000000

int main() {
int N, D;
double x = 0.99998000179983804;

double s, t;
s = t = 1.0;

while (s <= MAX_VAL && t <= MAX_VAL) {
if (fabs(x - s/t) < fabs(x - 1.*N/D)) {
N = s;
D = t;
}
if (s / t > x)
++t;
else
++s;
}

printf("%.18lf = %d/%d\n", x, N, D);
}

見事、元の分数の分子、分母を再現できる。

厳密な証明ではないが、数学的帰納法的なイメージで、上記の方法で適切なN, Dが求まることが分かる。(ざっくりした証明的なものです。。数学得意な人がいたら、厳密な証明方法教えてください。)

|x - N/D|が最小となるN, DをそれぞれN*, D*とおく。
i) (N* - 1, D*)、または、(N*, D* -1)が探索されれば、最適解にたどりつく。
ii) (N* - 1, D*) を (N~, D~) と改めて置きなおす。(N~-1, D~)、または、(N~, D~-1)が探索されれば、(N~, D~)に辿り着く。(N*, D*-1)に対しても同様。
iii)N, Dは1以上の整数なので、点列(D, D-1, D-2, ...)、および、(N, N-1, N-2, ..)は有限回の試行で、それぞれ1に収束する。
i)、ii)、iii)、より(N, D) = (1, 1)からはじめて、上記の方法を繰り返せば、最適解(N*, D*)に辿りつく。□

辿り着いた後も、N, D < Lの間は探索は続くので、最良解はメモリ上に保持しておかないといけないので注意。

2011年7月22日金曜日

ビットマップで画像処理

C++を使ってビットマップで画像処理をして遊んでみた。

実はC++の画像処理は昔仕事でちょっとだけやったことがあった。そのときは友達がつくったライブラリを使っていろいろ処理をしていたが、せっかくなので、自分でも作ってみた。

この神サイトが最強なので、ビットマップのデータ構造については割愛。

上記のサイトを参考にC++でビットマップ処理クラスを実装してみる。

まずは、構造体とかクラスの定義とか。

#include <cstdio>
#include <iostream>
#include <cstring>

#define FILE_HEADER_SIZE 14 // ファイルヘッダのサイズ
#define INFO_HEADER_SIZE 40 // 情報ヘッダのサイズ

using namespace std;

/*
* ファイルヘッダー構造体
*/
typedef struct FileHeader {
uint8_t data[FILE_HEADER_SIZE]; // 加工前データ(書き出しのため必要)
string fileType; // ファイルタイプ
int fileSize; // ファイルサイズ
} FileHeader;

/*
* 情報ヘッダー構造体
*/
typedef struct InfoHeader {
uint8_t data[INFO_HEADER_SIZE]; // 加工前データ(書き出しのため必要)
int infoHeaderSize; // 情報ヘッダのサイズ
int width; // 画像の幅
int height; // 画像の高さ
int clrPerPixel; // 1ピクセル当たりの色数
int dataSize; // 画像データのサイズ
} InfoHeader;

/*
* ピクセルの色情報構造体
*/
typedef struct Color {
int r; // 赤成分
int g; // 緑成分
int b; // 青成分
} Color;

int bit2Integer(uint8_t b1, uint8_t b2, uint8_t b3, uint8_t b4);

/*
* ビットマップ処理クラス
*/
class BitMapProcessor {
FILE *bmp; // ビットマップのファイルポインタ
uint8_t *img; // ビットマップデータ(加工用)
uint8_t *org; // ビットマップデータ(読み込み時)
FileHeader fHeader; // ファイルヘッダ
InfoHeader iHeader; // 情報ヘッダ

public:
BitMapProcessor() {
bmp = NULL;
img = NULL;
org = NULL;
};

~BitMapProcessor() {
fclose(bmp);
delete []img;
delete []org;
}
void loadData(string filename);
void dispBmpInfo();
void writeData(string filename);
Color getColor(int row, int col);
void setColor(int row, int col, int r, int g, int b);
void restore();
int height() { return iHeader.height; };
int width() { return iHeader.width; };
private:
void readFileHeader();
void readInfoHeader();
void readBmpData();

};

次に、定義したクラス、関数の実装。

/*
* 4ビット情報をInt整数値に変換
*/
int bit2Integer(uint8_t b1, uint8_t b2, uint8_t b3, uint8_t b4) {
return b1 +
b2 * 256 +
b3 * 256 * 256 +
b4 * 256 * 256 * 256;
}

/*
* ビットマップデータをロードする
*/
void BitMapProcessor::loadData(string filename) {
if (bmp != NULL)
fclose(bmp);

bmp = fopen(filename.c_str(), "rb");
if (bmp == NULL)
printf("ファイルオープンに失敗しました。\n");

readFileHeader();
readInfoHeader();
readBmpData();
}

/*
* ファイルヘッダを読む
*/
void BitMapProcessor::readFileHeader() {
uint8_t data[FILE_HEADER_SIZE];
size_t size = fread(data, sizeof(uint8_t), FILE_HEADER_SIZE, bmp);

memcpy(fHeader.data, data, sizeof(data));
fHeader.fileType = "";
fHeader.fileType += data[0];
fHeader.fileType += data[1];
fHeader.fileSize = bit2Integer(data[2], data[3], data[4], data[5]);
}

/*
* 情報ヘッダを読む
*/
void BitMapProcessor::readInfoHeader() {
uint8_t data[INFO_HEADER_SIZE];
size_t size = fread(data, sizeof(uint8_t), INFO_HEADER_SIZE, bmp);

memcpy(iHeader.data, data, sizeof(data));
iHeader.infoHeaderSize = bit2Integer(data[0], data[1], data[2], data[3]);
iHeader.width = bit2Integer(data[4], data[5], data[6], data[7]);
iHeader.height = bit2Integer(data[8], data[9], data[10], data[11]);
iHeader.clrPerPixel = bit2Integer(data[14], data[15], 0, 0);
iHeader.dataSize = bit2Integer(data[20], data[21], data[22], data[23]);
}

/*
* 画像データを読む
*/
void BitMapProcessor::readBmpData() {
if (img != NULL)
delete []img;

int sz = iHeader.dataSize;
img = new uint8_t [sz];
size_t size = fread(img, sizeof(uint8_t), sz, bmp);
if (size != sz)
printf("画像データ読み込みのサイズが矛盾しています。");

// バックアップ用にorgに画像データをコピー
if (org != NULL)
delete []org;
org = new uint8_t [sz];
memcpy(org, img, sz);
}

/*
* ビットマップ情報の表示(デバッグ用)
*/
void BitMapProcessor::dispBmpInfo() {
cout << "■ファイルヘッダ情報" << endl;
cout << "ファイルタイプ: " << fHeader.fileType << endl;
cout << "ファイルサイズ: " << fHeader.fileSize << endl;

cout << "■情報ヘッダ情報" << endl;
cout << "情報ヘッダサイズ: " << iHeader.infoHeaderSize << endl;
cout << "画像幅: " << iHeader.width << endl;
cout << "画像高: " << iHeader.height << endl;
cout << "1ピクセルあたりの色数: " << iHeader.clrPerPixel << endl;
cout << "画像データのサイズ: " << iHeader.dataSize << endl;
}

/*
* ビットマップデータのファイル書き出し
*/
void BitMapProcessor::writeData(string filename) {
FILE *out = fopen(filename.c_str(), "wb");

if (out == NULL)
printf("書き出し先のファイルを開けません。\n");

fwrite(fHeader.data, sizeof(uint8_t), FILE_HEADER_SIZE, out);
fwrite(iHeader.data, sizeof(uint8_t), INFO_HEADER_SIZE, out);
fwrite(img, sizeof(uint8_t), iHeader.dataSize, out);

fclose(out);
}

/*
* 指定されたピクセルの色を取得
*/
Color BitMapProcessor::getColor(int row, int col) {
if (row < 0 || row >= iHeader.height)
printf("getColor(): rowが範囲外です。\n");
if (col < 0 || col >= iHeader.width )
printf("getColor(): colが範囲外です。\n");

int width = 3 * iHeader.width;
while (width % 4) // ビットマップの1列は4の倍数ビットからなる
++width;

int bPos = row * width + 3 * col;
int gPos = bPos + 1;
int rPos = bPos + 2;

Color color;
color.r = img[rPos];
color.g = img[gPos];
color.b = img[bPos];

return color;
}

/*
* 指定されたピクセルに色を設定
*/
void BitMapProcessor::setColor(int row, int col, int r, int g, int b) {
if (row < 0 || row >= iHeader.height)
printf("getColor(): rowが範囲外です。\n");
if (col < 0 || col >= iHeader.width )
printf("getColor(): colが範囲外です。\n");

int width = 3 * iHeader.width;
while (width % 4) // ビットマップの1列は4の倍数ビットからなる
++width;

int bPos = row * width + 3 * col;
int gPos = bPos + 1;
int rPos = bPos + 2;

img[rPos] = r;
img[gPos] = g;
img[bPos] = b;
}

/*
* ビットマップデータを加工前に復元する
*/
void BitMapProcessor::restore() {
memcpy(img, org, iHeader.dataSize);
}

あとは、適当にテストして遊ぶ。

/*
* テスト用関数(1)モノクロ化
*/
void twoTone(BitMapProcessor *bmp) {
for (int i = 0; i < bmp->height(); i++)
for (int j = 0; j < bmp->width(); j++) {
int ave = 0;
ave += bmp->getColor(i, j).r;
ave += bmp->getColor(i, j).g;
ave += bmp->getColor(i, j).b;
ave /= 3;

bmp->setColor(i, j, ave, ave, ave);
}
}

/*
* テスト関数(2)指定範囲の切り取り
*/
void extractArea(BitMapProcessor *bmp, int r0, int r1, int c0, int c1) {
for (int i = 0; i < bmp->height(); i++)
for (int j = 0; j < bmp->width(); j++) {
if (r0 <= i && i <= r1 && c0 <= j && j <= c1)
continue;
bmp->setColor(i,j, 255, 255, 255);
}
}

/*
* テスト関数(3) 色の反転
*/
void invert(BitMapProcessor *bmp) {
for (int i = 0; i < bmp->height(); i++)
for (int j = 0; j < bmp->width(); j++) {
int ave = 0;
int r = bmp->getColor(i, j).r;
int g = bmp->getColor(i, j).g;
int b = bmp->getColor(i, j).b;

bmp->setColor(i, j, 255-r, 255-g, 255-b);
}
}

/*
* テスト関数(4)モザイク化
*/
void mosaic(BitMapProcessor *bmp, int level) {
if (level <= 0)
level = 1;

for (int i = 0; i < bmp->height(); i+=2*level)
for (int j = 0; j < bmp->width(); j+=2*level) {
int r = 0;
int g = 0;
int b = 0;
int cnt = 0;

for (int x = -level; x <= level; x++)
for (int y = -level; y <= level; y++) {
int xt = i + x;
int yt = j + y;

if (xt < 0 || yt < 0 || xt >= bmp->height() || yt >= bmp->width())
continue;
++cnt;
r += bmp->getColor(xt, yt).r;
g += bmp->getColor(xt, yt).g;
b += bmp->getColor(xt, yt).b;
}

r /= cnt;
g /= cnt;
b /= cnt;

for (int x = -level; x <= level; x++)
for (int y = -level; y <= level; y++) {
int xt = i + x;
int yt = j + y;

if (xt < 0 || yt < 0 || xt >= bmp->height() || yt >= bmp->width())
continue;

bmp->setColor(xt, yt, r, g, b);
}
}
}

/*
* メイン処理
*/
int main() {
BitMapProcessor bmp;

// ビットマップデータのロード
bmp.loadData("kii.bmp");

// ビットマップ情報の表示
bmp.dispBmpInfo();

// テスト1. モノクロ化
twoTone(&bmp);
bmp.writeData("kii_1.bmp");
bmp.restore();

// テスト2. 指定範囲の切り出し
extractArea(&bmp, 200, 300, 100, 180);
bmp.writeData("kii_2.bmp");
bmp.restore();

// テスト3. 色反転
invert(&bmp);
bmp.writeData("kii_3.bmp");
bmp.restore();

// テスト4. モザイク
mosaic(&bmp, 30);
bmp.writeData("kii_4.bmp");
bmp.restore();

return 0;
}


こんな感じ。これはなかなか面白い!

かわいいっ!


やっぱ、かわいい!

ちょっとこわいけど、かわいい!

最高!
それでも、かわいい!

2011年7月18日月曜日

Haskell勉強記(8)

前回分からなかったfoldlをfoldrで書きかえる方法について。
duality theorems concerning foldl and foldrの一つらしい。
main = print $ foldl (++) [] ["1", "2", "3", "4", "5"]
は、foldrを用いて以下のように書ける。
main = print $ foldr (\x g n -> g (n ++ x)) id ["1", "2", "3", "4", "5"] []
関数適用がもっとも結合度が高いことと、関数適用はleft-associativeであることを考えると、上式は以下と等価。
main = print $ foldr (\x g n -> g (n ++ x)) id ["1", "2", "3", "4", "5"] $ []
ちょっとだけ分かりやすく書くと、
func = foldr (\x g n -> g (n ++ x)) id ["1", "2", "3", "4", "5"]
main = print $ func []
なるほど、foldrは、値と関数を受け取って新しい関数を返す関数であることが分かる。やってるのは、こんなこと。
func = (\x g n -> g (n++x)) [1] $
(\x g n -> g (n++x)) [2] $
(\x g n -> g (n++x)) [3] $
(\x g n -> g (n++x)) [4] $
(\x g n -> g (n++x)) [5] $ id

main = print $ func []

それ自体を返す関数を受け取って、5を末尾に付ける関数を作る。それを受け取って、4を末尾に付けた後に、5を末尾に付ける関数を作る。さらにそれを受け取って、・・・・・という具合。

一般に以下の書き換えが可能。


オイラー関数の公式の確率的意味

最近、オイラー関数(Euler's totient function)について少し復習してみた。新しい発見があったのでメモしておく。

オイラー関数とは、



と素因数分解できるある自然数nに対して、


と表され、これはn未満の自然数でnと互いに素なものの数を表しています。

 上記の公式の証明は、以下3つの補題から簡単にできます。







1番目の補題は、pが素数なことから自明。2番目の補題は、p^k未満の数でpを素因数にもつものを数えれば導出できます。3番目の補題は、aと因数を共有しないもの、bと因数を共有しないものを、それぞれaとbの法を元に列挙して、中国の剰余定理を使えば証明できます。
これら3つの補題からオイラー関数の公式は導出可能ですが、詳細は省略。

面白いと思ったのは、オイラー関数を
(n未満のnと互いに素な自然数の個数)
= (1-nまでの自然数の個数)
* (1-nまで自然数のp1で割れない確率)
* (1-nまで自然数のp2で割れない確率)
.....
* (1-nまで自然数のpkで割れない確率)

と捉えることができるということです。確かに言われてみればその通りですが、この発想は面白いと思いました。

2011年7月14日木曜日

C++でハッシュのマップを実装する

C++のset/mapは、二分木で実装されている。要素数が大きくなると重くなってしまうので(平衡木の回転の処理が重たい?)、ハッシュを使ってmapみたいなものを実装してみた(挿入、参照のみ)。

template <class K, class V>
class hashMap {
static const unsigned int DEFAULT_SIZE = 1000000;
vector<pair<K, V> > *contents;
unsigned int _size;
unsigned int (*_hashValue) (K, unsigned int);

public:
hashMap(unsigned int (*func) (K, unsigned int) , unsigned int size = DEFAULT_SIZE) {
_hashValue = func;
_size = size;
contents = new vector<pair<K, V> > [_size];
}

~hashMap() {
delete [] contents;
}

bool hashElement(K s) {
unsigned int h = _hashValue(s, _size);

for (int i = 0; i < (int)contents[h].size(); i++) {
if (contents[h][i].first == s)
return true;
}

return false;
}

void insert(K x, V n) {
if (hashElement(x))
return;

unsigned int h = _hashValue(x, _size);
contents[h].push_back(make_pair(x, n));
}

pair<K, V> *getElement(K x) {
unsigned int h = _hashValue(x, _size);
for (int i = 0; i < (int)contents[h].size(); i++)
if (contents[h][i].first == x)
return &contents[h][i];

return NULL;
}

void printHashBalance() {
int cnt = 0;
double ret = 0.0;
int worst = 0;

for (int i = 0; i < (int)_size; i++) {
if (contents[i].size()) {
++cnt;
ret += contents[i].size();
worst = max<int>(worst, contents[i].size());
}
}

if (!cnt)
cerr << "Hash is empty" << endl;
else {
cerr << "collisions occurred(average): " << (ret/cnt) << endl;
cerr << "collisions occurred(worst): " << worst << endl;
}
}
};


これで、参照/挿入はO(1)でできるはず。衝突がどれくらい起きているかを調べるためにprintHashBalance()という関数を定義した。コンストラクタには、使用するハッシュ関数のポインタと、ハッシュ値のサイズを渡す。

適当にテストコードを書いて、標準関数のmapと比較する。WORDS_NUM個のランダムな文字列を、それぞれのコンテナにぶちこんで、最頻出の文字列を求める。

#define WORDS_NUM 10000

/*
* generate random strings of length between 3 and 8 inclusive
*/
vector<string> randomGen() {
vector<string> ret;

for (int i = 0; i < WORDS_NUM; i++) {
int len = 3 + rand() % 6;

char str[len+1];
for (int j = 0; j < len; j++)
str[j] = 97 + rand() % 26;
str[len] = '\0';

ret.push_back(str);
}

return ret;
}
/*
* timer
*/
class Timer {
double start_t;
public:
void start() {
cout << "Timer Start." << endl;
start_t = gettimeofday_sec();
}

void end() {
cout << "Timer stopped." << endl;
cout << "Time Elapsed: " << gettimeofday_sec() - start_t << endl;
}
};

//#define _MAP
#define _HASH

unsigned int stringHash(string str, unsigned int sz) {
unsigned int ret = 0;

for (int i = 0; i < (int)str.size(); i++)
ret = 26 * ret + str[i] - 97;

return ret % sz;
}

int main() {
vector<string> words = randomGen();
Timer timer;
timer.start();

#ifdef _MAP
map<string, int> wc;
int cnt = 0;
string best;

for (int i = 0; i < (int)words.size(); i++) {
if (wc.count(words[i])) {
int n = ++wc[words[i]];

if (n > cnt) {
cnt = n;
best = words[i];
}
}
else
wc[words[i]] = 1;
}
cout << cnt << ": " << best << endl;
#endif

#ifdef _HASH
hashMap<string, int> wc(stringHash, WORDS_NUM);
int cnt = 0;
string best;

for (int i = 0; i < (int)words.size(); i++) {
if (wc.hashElement(words[i])) {
int n = ++wc.getElement(words[i])->second;

if (n > cnt) {
cnt = n;
best = words[i];
}
}
else
wc.insert(words[i], 1);
}
cout << cnt << ":" << best << endl;
wc.printHashBalance();
#endif

timer.end();

return 0;
}

以下比較結果。HashMapの_sizeは入力文字数と同じサイズにした。一つのハッシュ値に格納される入力値の種類は平均1.45くらいだった。

WORDS_NUM100,0001,000,00010,000,000
Map [s]0.3464.99966.71
HashMap [s]0.22.05220.274

もっと速くなると思ったが、どこがネックなんだろう。オーダー的にはハッシュ値の計算が重いけど、実際はvectorのメモリ関連処理が遅そうな気がする。プロファイルとってみるかー。

2011年7月13日水曜日

Haskellで有名アルゴリズム実装(5)

最近英語力が落ちてる感があるので、英語で書く。

The fifth trial of implementing a classic algorithm with Haskell!
Today I'll try Pascal' Triangle.
pascalTriangle x = x : pascalTriangle y
where y = zipWith (+) (0:x) (x++[0])
Nothing special. Thanks to lazy evaluation, you can define an infinite list of Pascal's Triangle.

Here's a similar problem:
Let's consider how many times it lands heads if you toss a coin. More specifically, what the probability of the coin landing heads k times out of n tosses?

It seems it's nothing to do with Pascal's Triangle. But it actually is.
prbTree x = x : prbTree y
where y = zipWith (+) (map (*0.5) (0:x)) (map (*0.5) (x++[0]))
The tree structure is just like Pascal's Triangle.
In the above source, I used floating point numbers. So they don't sum up to 1.0 due to the limitation of accuracy of the type.

You can also use rational numbers, which is in Data.Ratio module, instead of floating point numbers. In that case, they sum up to 1.0 when added together.
import Data.Ratio

half = (*) (1%2)
prbTree x = x : prbTree y
where y = zipWith (+) (map half $ 0:x) (map half $ x++[0])

2011年7月10日日曜日

クイックソート VS マージソート

今日は、軽い実験。Haskellの練習がてら、quick sortとmerge sortを比較する。

まず、quick sort。
-- quick sort
qsort :: (Ord a) => [a] -> [a]
qsort [] = []
qsort (x:xs) = let left = filter (<x) xs
right = filter (>=x) xs
in qsort left ++ [x] ++ qsort right
Haskellは美しい。次に、merge sort。
msort :: (Ord a) => [a] -> [a]
msort [x] = [x]
msort [x, y] = [min x y, max x y]
msort xs = let len = length xs `div` 2
pre = msort $ take len xs
suf = msort $ drop len xs
in merge pre suf

merge [] [] = []
merge x [] = x
merge [] y = y
merge (x:xs) (y:ys)
| x <= y = x : merge xs (y:ys)
| otherwise = y : merge (x:xs) ys
もうちょっと綺麗に書けそうだが、こんなもんだろう。
適当に乱数をソートしてみる。
getRandoms :: Int -> [Int]
getRandoms n = take n $ randomRs (1, 100000) $ mkStdGen 0
n=100,000で比較。
quick sort 0.653[s]
merge sort 0.891[s]

quick sortの弱点を突いてみる。ソート済みのリストを入力として与えてみる。[1..100000]で比較。
quick sort 815.877[s]
marge sort 0.465[s]

quick sortはpivot選択すれば改良はできるけど、それでも理論上の最悪計算量はO(n^2)。merge sortは常に平衡二分木を形成するように分割するので、計算量は常にO(n log n)。ランダムなデータに対してはquick sortが速い場合が多いようだが、実用上はmerge sortを使用した方が安心な気がする。

2011年7月7日木曜日

モナドを実装する

 今日はモナド実装にチャレンジ。
モナドの機能は、contextualな値に関数を適用すること。
functorやapplicativeも似たようなことを行うが、モナドのいいところは演算子のつなぎ目をきれいに書けるところ。

 とりあえず、書いてみる。まずモナドクラスのインスタンスにする型を定義する。
Maybe型を応用して、即成功、即失敗、判断を後回しにする、という3つのcontextを持つ型を定義する。
data Check a = YES | NO | Unknown a deriving (Show, Eq)
 この型を利用するために以下の関数を定義する。
passTest f x  = if f x == True then YES else Unknown x
failTest f x = if f x == True then NO else Unknown x
determine f x = if f x == True then YES else NO
 さらに、この型をモナドインスタンスとして実装する。
instance Monad Check where
return x = Unknown x
YES >>= _ = YES
NO >>= _ = NO
Unknown x >>= f = f x
 これを使うと、ややこしい制御フローがきれいに書ける。例えば、うるう年判定をするプログラム。やたらとif-elseが多い。
isLeapYear x = if x `mod` 400 == 0
then True
else if x `mod` 100 == 0
then False
else if x `mod` 4 == 0
then True
else False
 モナドを使うときれいに書ける。
divisable y x = x `mod` y == 0
isLeapYear x = return x >>= passTest (divisable 400)
>>= failTest (divisable 100)
>>= determine (divisable 4)
 モナドは、制御フローを表に出さずに演算を繋げることができるデザインパターンと考えることもできる。見てのとおり制御フローをcomposableな要素として部品化することができ、if-else-break文のボイラープレート除去を行うことができる。条件が増えたら、新しいバインドで関数をつないであげれば、ifとかelseとかbreakとか書かなくてもよい。(composableうんぬんは、tanakhさんがつぶやいていた話です。多分もっと高尚なお話だったと思いますが、自分が理解した限りだとこんな話。。のはず。)

functorやapplicativeと比べて、モナドがやたらと注目されている理由は、この3つの中で最も抽象度が高いからという理由もあるが、上のようなデザインパターンとして使用できるからというのが大きいのだと思う。

最後に、今回実装したモナドがモナド則を満たしていることを確認する。
適当に値を与えて確認する。1つ1つやってもいいけど、せっかくなので最近覚えたlist applicativeを使って複数個のパターンを一気に確認する。
nums  = [-5, 0, 5]
ck = [YES, NO, Unknown (-5), Unknown 0, Unknown 5]
funcs = [passTest (divisable 400),
failTest (divisable 100),
determine (divisable 4)
]

assertLeftIdentity = (\f x -> (return x >>= f) == f x) <$> funcs <*> nums
assertRightIdentity = (\m -> (m >>= return) == m) <$> ck
assertAssociativity = (\f g m -> ((m >>= f) >>= g)
== (m >>= (\x -> f x >>= g))) <$> funcs <*> funcs <*> ck

イケてた。モナドイケイケ~~。

2011年7月4日月曜日

Haskellで有名アルゴリズム実装(4)

第4弾。今日は、素数列挙のアルゴリズム。
ある数が素数かどうかの判定は、擬多項式時間でできます。ある区間[1, n]内の素数を列挙するには、普通にやるとO(n*sqrt(n))かかります。
エラトステネスの篩というアルゴリズムを使用すると、O(n)で計算できます。

では、実装。
Haskellのエラトステネスは美しいです。
eratos []     = []
eratos (x:xs) = x : eratos [y | y <- xs, mod y x /= 0]
これで十分な場合がほとんどですが、100万までの素数列挙とかすると、上の実装ではたえられません。エラトステネスが速いのは、ある素数pの倍数を2*p, 3*p,..とpおきに消していけるからで、上の方法だとすべての要素をなめているので遅いです。加えてリストは参照にO(n)かかるので遅いです。

ある素数に対して必要な要素だけ参照する+高速に参照することができるデータ構造に変換する。というのが必要です。100万くらいの素数列挙なら、Setを使えば十分でしょう。あとは、thunkがメモリを圧迫するのを防ぐ+遅延評価によるスタックオーバーを防ぐため正格評価を用います。
primes :: [Integer] -> [Integer]
primes = toList . (seive 2) . fromList

seive :: Integer -> Set Integer -> Set Integer
seive p xs
| p*p > findMax xs = xs
| otherwise = seive (p+1) $ foldl' (\acc x -> delete x acc) xs [2*p,3*p..q]
where q = findMax xs
ただ、上の実装でも本来のO(n)の速度には届かないので、より大きな範囲([1, 10000000]など)の素数列挙をしたい場合は、工夫が必要です。Setは参照がO(log n)なので、ランダムアクセスが可能なデータ構造を利用する + なるべく遅延評価の機構をいれないようにするで達成できるかと思います。

2011年7月2日土曜日

Haskellでメモ化する(中級)

 前回のpostで、Haskellでメモ化する方法について書きました。そのとき、疑問点として、「なんで、メモ化と不動点が繋がってるねん!なんでやねん!」というのがありました。その疑問が解決したので中級編を書きます。

 どうやら不動点コンビネータを利用してやりたいことは、「任意の関数が与えられたときに、その計算をメモ化する関数を作成する」ということのようです。この発想にはびっくりです。どれだけ抽象度の高いことをやるんだよ。という感じです。imparative、あんど、object-orientedな思想に育てられた(毒された?)普通のプログラマーである私には思いもよらない発想でした。。

 まずは、メモ化と、不動点コンビネータについて少しおさらいをします。フィボナッチ数列をメモ化再帰で解くプログラムを書いてみます。
memoFib =
let fib 0 = 1
fib 1 = 1
fib n = memoFib (n-1) + memoFib (n-2)
in (map fib [0..] !!)
これは、list-comprehensionを用いて以下のように書くこともできます。(こっちの方が分かりやすいかも??)
memoFib = [fib x | x <- [0..]]
where fib 0 = 1
fib 1 = 1
fib n = memoFib !! (n-1) + memoFib !! (n-2)
と、とりあえず、関数がgivenであれば、その関数をメモ化してあげることは出来ます。

 次に、不動点コンビネータについて。Haskellでは、Control.Monad.Fixに不動点コンビネータが定義されています。これを使ってもいいですが、勉強のために自分で実装してみます。再帰関数を不動点コンビネータに書きかえる方法についてはここのページがとても分かりやすいです。
フィボナッチ関数を不動点コンビネータで書きなおします。フィボナッチ数列は、ラムダ式を用いて
fib = (\x -> if x <= 1 then 1 else fib (x-1) + fib (x-2))
と書けます。さらに、再帰で呼び出される関数を引数で与えることにすると、
fib = (\f x -> if x <= 1 then 1 else f (x-1) + f (x-2)) fib
と書けます。上の式をみると、フィボナッチ関数は
(\f x -> if x <= 1 then 1 else f (x-1) + f (x-2))
の不動点となっています。つまりフィボナッチ関数は、
fix (\f x -> if x <= 1 then 1 else f (x-1) + f (x-2))
と書けるというわけです。(詳しくは上のお勧めページ参照。)

 と、準備はここまで。いよいよ任意の関数をメモ化関数に書きかえる関数を定義するわけですが、任意の関数というのは、Haskellビギナーの私にはちと難しいので、(Int -> Int)型の任意の関数をメモ化することにします。

 まず、引数にそのままメモ化したい再帰関数を渡してしまうと意味がない(メモ化機構を入れられない)ので、メモ化したい関数を不動点にとるような関数を引数として渡してあげます。
fix f = f (fix f)
とすると、fを渡すという意味です。
ここまで来ると、なんとなくできそうな感じがします。なんとなく実装したのが、これ。
memorizeCalc f = ([f (memorizeCalc f) x | x <- [0..]] !!)
計算自体は正しいですが、メモ化が機能していないようです。前にも書きましたが、Haskellの場合、グラフ簡約のメモ化(ここでいうメモ化は、グラフ簡約されて同じ計算がされないという意味のメモ化)が働く場合と働かない場合があります。この形は前に見たことがあります。おそらく、引数をやめて実体を無限リストとして定義してあげればグラフ簡約が働くはず。で、書いたのがこちら。
memorizeCalc f = memo
where memo = ([f memo x | x <- [0..]] !!)
これで、「与えられた関数をメモ化関数に変換する関数」が出来ました。実際に確かめてみます。
fib' = (\f x -> if x <= 1 then 1 else f (x-1) + f (x-2))
trippleFib' = (\f x -> if x <= 2 then 1 else f (x-1) + f (x-2) + f (x-3))

main = do print $ memorizeCalc fib' 300
print $ memorizeCalc trippleFib' 300

値も正しく、高速化もされていることが確認できました。

 次のテーマは、
  • (Int -> Int)型以外の関数のメモ化
  • 高速化
です。高速化については、無限平衡二分木を用いた方法が紹介されています。また、Stateモナドを用いた方法があるようです。他にもData.ArrayのArrayが使えるんじゃないかなあ、とにらんでいます。。この辺が出来たら、上級編を書きたいです。