(相補)誤差関数と標準正規分布の関係

math

まぁどこにでも書いてあるような話ですが……。
標準正規分布の分布関数 \[ F(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} \exp(-\frac{t^2}{2}) dt \]

が欲しかったのだがC++の標準ライブラリには入っていない。しかし<cmath>にある相補誤差関数erfcを利用してこれを計算することができる。
誤差関数\({\rm erf}(x) \)は \[ {\rm erf}(x) \equiv \frac{2}{\sqrt{\pi}} \int_0^x \exp(-t^2) dt \]

と定義される。相補誤差関数は\({\rm erfc}(x) = 1 - {\rm erf}(x) \)であり、 \begin{eqnarray} {\rm erfc}(x) &\equiv& 1 - {\rm erf}(x) \\ &=& \frac{2}{\sqrt{\pi}} \int_0^\infty \exp(-t^2) dt - \frac{2}{\sqrt{\pi}} \int_0^x \exp(-t^2) dt \\ &=& \frac{2}{\sqrt{\pi}} \int_x^\infty\ \exp(-t^2) dt \\ &=& \frac{2}{\sqrt{\pi}} \int_{-\infty}^{-x}\ \exp(-t^2) dt. \end{eqnarray}

途中、ガウス積分
\[ \int_{-\infty}^\infty \exp(-ax^2) dx = \sqrt{\frac{\pi}{a}} \]

および\(\exp(-t^2) \)の対称性を用いた。ここで\( t = \frac{\xi}{\sqrt{2}} \)とおけば

\[ {\rm erfc}(x) = \frac{2}{\sqrt{\pi}} \int_{-\infty}^{-\sqrt{2}x}\ \exp(-\frac{\xi^2}{2}) d(\frac{\xi}{\sqrt{2}}). \]

したがって分布関数は \begin{eqnarray} F(x) &=& \frac{1}{2} {\rm erfc}(-\frac{x}{\sqrt{2}}) \\ &=& \frac{1}{2}(1 - {\rm erf}(-\frac{x}{\sqrt{2}})) \\ &=& \frac{1}{2}(1 + {\rm erf}(\frac{x}{\sqrt{2}})) \end{eqnarray}

と表現できる。よって

  double normal (double d) {
    return ( 0.5 * erfc ( -d * M_SQRT1_2 ) );
  }

などとすればよいのでした。

参考文献

debian 7.8 の vim で :filetype が通らない話

vim debian

Windows でコーディングするのがあまりにも辛い (MinGW とか Cygwin 的な意味で) ので debian 7.8 をインストールした．で，github から dotfiles を持ってきて :NeoBundleInstall しようとしたら通らないのでエラーコードを見たらそもそも .vimrc の頭にある :filetype off が通ってない．どうやら debian のデフォルトの vim は vim-tiny というやつでいろいろなオプションが無効になっているらしい．バージョンもちょっと古いので最新のソースを持ってきてビルドした．

% hg clone https://vim.googlecode.com/hg/ vim
% cd vim
% ./configure --with-features=huge #後ろにお好みのオプション
% make
% sudo make install

これでとりあえず無事に動いた．

参考

Vim の種類 (Vim family) - Qiita

vim-jp » Linuxでのビルド方法

Eratosthenes の篩

cpp

前回のエントリの続きでまた素数の話．
Team-lablog » 1000万個目の素数を超高速に出力せよ
これに触発されて 1 000 万番目までの素数の出力とその高速化に挑戦してみる．
まずは前回のソースコードをちょっと手直しして，

#include <iostream>
#include <ctime> //実行時間の測定
using namespace std;
#define MAX 180000000

int main(){
	int count = 0;
	clock_t begin, finish;
	begin = clock();
	bool *num;
	num = new bool[MAX];
	for(int i = 0; i < MAX; i++){
		num[i] = true;
	}

	num[0] = false;
	num[1] = false;
	
	for(int i = 2; i * i <= MAX; i++){
		for(int j = i * 2; j <= MAX; j += i){
			if(!num[j]){
				continue;
			}
			num[j] = false;
		}
	}

	for(int i = 0; i < MAX; i++){
		if(num[i]){
			cout << i << " ";
			count++;
			if(count == 10000000){
				cout << endl << "Finished!" << endl;
				break;
			}
		}
	}

	finish = clock();
	cout << (finish - begin) / 1000 << "ms" << endl;

	delete[] num;
	return 0;
}

1 000 万番目の素数は 179 424 673 らしい*1ので最大値は 180 000 000 とした．
プロセッサは Core i5 の 1.3GHz．
で，実行結果．
f:id:kailily:20140809054132p:plain
何も工夫してないのでこんなものか．
偶数を飛ばすように書き換える．

#include <iostream>
#include <ctime>
using namespace std;
#define MAX 180000000

int main(){
	int count = 0;
	clock_t begin, finish;
	begin = clock();
	bool *num;
	num = new bool[MAX];

	for(int i = 2; i < MAX; i += 2){
		num[i] = false;
		num[i + 1] = true;
	}

	num[0] = false;
	num[1] = false;
	num[2] = true;

	for(int i = 3; i * i <= MAX; i += 2){
		for(int j = i * i; j <= MAX; j += i * 2){
			if(!num[j]){
				continue;
			}
			num[j] = false;
		}
	}

	for(int i = 0; i < MAX; i++){
		if(num[i]){
			cout << i << " ";
			count++;
			if(count == 10000000){
				cout << endl << "Finished!" << endl;
				break;
			}
		}
	}

	finish = clock();
	cout << (finish - begin) / 1000 << "ms" << endl;

	delete[] num;
	return 0;
}

実行結果．
f:id:kailily:20140809134642p:plain
かなり速くなった．
まだ無駄な繰り返しがあるので速く出来るとは思うが，思い付かなかったので良さそうな方法を見つけたらまた書く．

*1:http://primes.utm.edu/curios/page.php/179424673.html

AOJ 1004

cpp aoj

簡単かと思ったけど意外と手こずってしまった．
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1004
入力 n に対して，1 から n までの列とそれをひっくり返した列 (n から 1 まで) を重ねた上下をペアとして，ふたつとも素数であるようなペアの数を数える．
最初は素数を判定させる関数を作っていちいちそれに放り込む作戦をとったが，時間切れになってしまうので boolean 型の配列そのまま使った．

#include <iostream>
using namespace std;

void prime();
bool num[10001];

int main(){
	int m, count;
	prime();
	while(cin >> m){
		count = 0;
		for(int i = 1; i <= m; i++){
			if(num[i] && num[m - i + 1]){
				count++;
			}
		}
		cout << count << endl;
	}
	return 0;
}

void prime(){
	for(int i = 0; i < 10001; i++){
		num[i] = true;
    	}
    	num[0] = false;
    	num[1] = false;
	for(int i = 2; i <= 100; i++){
        	for(int j = i*2; j <= 10001; j += i){
            		if(!num[j]){
                		continue;
            		}
            		num[j] = false;
        	}
    	}
}

素数を判定するプログラムはたまに必要になるのでメモしておく．
とりあえず 1 000 000 以下の素数について (AOJ 0009 の使い回しなので)．

#include <iostream>
using namespace std;

int main(){
	int n;
	bool num[1000001];
	for(int i = 0; i < 1000001; i++){
		num[i] = true;
	}
	num[0] = false;
	num[1] = false;
	for(int i = 2; i <= 1000; i++){
		for(int j = i * 2; j <= 1000000; j += i){
			if(!num[j]){
				continue;
			}
			num[j] = false;
		}
	}
	while(cin >> n){
		if(num[n]){
			cout << n << " is prime." << endl;
		}
		else{
			cout << n << " is not prime." << endl;
		}
	}
	return 0;
}