## 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月30日土曜日

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

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

「(省略) これは、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;}

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:
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日火曜日

### 無理数を有理数で近似

N/D > x なら 探索点を N/(D+1)に移動、そうでなければ探索点を(N+1)/Dに移動。これを逐次繰り返していき、最良点を保持すれば答えが出る。

#define MAX_VAL 1000000int 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);}

|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*)に辿りつく。□

## 2011年7月23日土曜日

### ビットマップで画像処理

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; fHeader.fileType += data; fHeader.fileSize = bit2Integer(data, data, data, data);}/** 情報ヘッダを読む*/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, data, data, data); iHeader.width          = bit2Integer(data, data, data, data); iHeader.height         = bit2Integer(data, data, data, data); iHeader.clrPerPixel    = bit2Integer(data, data, 0, 0); iHeader.dataSize       = bit2Integer(data, data, data, data);}/** 画像データを読む*/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 << "１ピクセルあたりの色数: " << 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月19日火曜日

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"] []

main = print $foldr (\x g n -> g (n ++ x)) id ["1", "2", "3", "4", "5"]$ []
ちょっとだけ分かりやすく書くと、

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

## 2011年7月5日火曜日

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

では、実装。
eratos []     = []eratos (x:xs) = x : eratos [y | y <- xs, mod y x /= 0]