PythonとCでモンテカルロ <訂正版>

以前の記事(PythonとCでモンテカルロ法を用いて円周率を求めてみた)に重大なミスがあることが発覚しました。
結論から言うと、pythonとC言語で60倍近くも速度に差があると書きました。
あれは嘘です。
今回は、前回のプログラムに対して訂正を加えてより正確に”計算速度の差“を比べてみました。

 

前回プログラムの問題点

前回記事を書いた時にうすうす気づいていたのですが…C言語とPythonで乱数を生成する際に大きな差があるようです。詳しくは参考記事を読んでいただきたいのですが、ちょっと説明をすると、C言語の乱数は”線形合同法“という比較的簡単な方法で生成されています。対してPythonは”メルセンヌツイスタ“と呼ばれる複雑な計算をしてより精度の高い(乱数に近い?)乱数を生成しています。その乱数を生成する仮定で実行結果に大きく影響を与えていたことが発覚しました。

前回の計算結果を見てみるとわかるのですが、同じ回数乱数を生成して計算を行っているのに円周率はPythonの方がより正確です。

(1億回ループ時)
C : 3.145632
Python : Pi:3.141701  <- C言語の計算結果よりも円周率に近い
(正解 : 3.14159265…)

これはどうやら偶然じゃなく、乱数の精度の影響だったようです。。。

 

今回のプログラム

今回のプログラムでは、乱数生成時に不公平が内容に乱数をファイルから読み込んで、それに対して計算を行いました。
さらに、プログラムの実際に計算の部分について比較をしたかったので、プログラム中で時間を計測しました(ファイルの読み込み時は含まない)

ソースコード

C

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define RAND_SIZE 100000000
#define FILE_NAME "rand.txt"

typedef struct{
 float x;
 float y;
}pos;

int chk_pos( double, double, int);
void data_set( pos*, char*, int);

int main(void){
 int cnt, i;
 double start, loop_time;
// pos rand_data[RAND_SIZE];
// メモリを動的確保にしないと、メモリ確保でsegmentation faultで落ちる
 pos *rand_data2 = (pos *)malloc(sizeof(pos) * RAND_SIZE);

data_set( rand_data2, FILE_NAME, RAND_SIZE);

start = clock();
 cnt = 0;
 for(i=0; i<RAND_SIZE; i++){
  cnt += chk_pos( rand_data2[i].x, rand_data2[i].y, 1);
 }
 loop_time = clock() - start;

free(rand_data2);

printf("loop:%d pi:%f\n", RAND_SIZE, (double)cnt / RAND_SIZE * 4);
 printf("loop time : %f[s]\n", loop_time / CLOCKS_PER_SEC);

return 0;
}

// 中央からの距離が1以内であるか
int chk_pos( double x, double y, int r){
 int ret;
 double dist_r;

dist_r = sqrt(pow(x,2) + pow(y,2));

if (dist_r < r) ret = 1;
 else ret = 0;

return ret;
}

// ランダムデータの読み込み
void data_set(pos *data, char* file_name, int size){
 FILE *fp;
 char s[256];
 int i;
 float pos_x, pos_y;

if((fp = fopen(file_name, "r")) == NULL){
 printf("cannot open file\n");
 exit(EXIT_FAILURE);
 }

i=0;
 while(fgets( s, 256, fp) != NULL){
 if(i >= size) break;
 sscanf(s,"%f,%f", &pos_x, &pos_y);
 data[i].x = pos_x, data[i].y = pos_y;
 i++;
 }

fclose(fp);
}

 

Python

# ファイルから読み込み機能のついた円周率計算1
import random
import math
import time

# 定数
RAND_SIZE = 100000000
FILE_NAME = "rand.txt"

# ランダムのデータに関するクラス
class random_data:
 rand_data = []
 # インスタンスを作成した時点でデータを入力する
 def __init__( self, input_file_name, size):
 i = 0
 for line in open(input_file_name):
  if (i >= size):
   break
  line = line.replace("\n","")
  str_tmp = line.split(",")
  self.rand_data.append([float(str_tmp[0]),float(str_tmp[1])])
  i += 1

# データ表示
 def print_data( self):
 for data in self.rand_data:
  print("x : %f, y : %f" % ( data[0], data[1]))

# 中央からの距離が1以内であるか
def chk_pos(x,y,r):
 dist_r = math.sqrt( math.pow(x,2) + math.pow(y,2))
 if dist_r < r:
  ret = True
 else:
  ret = False
 return ret

data = random_data( FILE_NAME, RAND_SIZE)
# data.print_data();

cnt = 0
i = 0
# start = time.time();
start = time.clock();
while i < RAND_SIZE:
# x = random.random()
# y = random.random()
# if chk_pos(x,y,1):
 if chk_pos( data.rand_data[i][0], data.rand_data[i][1],1):
 cnt += 1
 i += 1
# loop_time = time.time() - start
loop_time = time.clock() - start

print("loop:%d Pi:%f" % ( RAND_SIZE, cnt / RAND_SIZE * 4))
print("loop time : %f[s]" % loop_time)

C言語のプログラムについて、メモリを動的確保しています。これは、コメントアウトでも書いてありますが、配列を宣言するさいにrand_data2[100000000]ってやると、それだけでsegmentation faultを出してプログラムが落ちたからです。
(ちなみに、rand_data2[1000000]までならとりあえず動作した)
謎ですねぇ。。。まぁ、とりあえずそーゆーことだったんで動的確保にしています。

プログラムについて、何かミスがあればTwitterかコメントで教えて頂けるとありがたいです。
(変数名や関数名にセンスを感じないって指摘はNGで)

 

結果

さて、お待ちかねの実行結果です。実行時間もついでに計測したので、timeをつけて実行しています。

# C言語
time ./pi2
--> loop:100000000 pi:3.141572
--> loop time : 1.476243[s]
--> ./pi2 97.26s user 3.97s system 96% cpu 1:44.72 total
# Pythn
time python pi2.py
--> loop:100000000 Pi:3.141572
--> loop time : 236.437197[s]
--> python pi2.py 635.43s user 950.89s system 45% cpu 57:35.71 total

わぉ。C言語が1.48秒に対して、Pythonが236.44秒…。
計算すると、なんと160倍C言語の方が速いって結果に！！！
プログラムの書き方で結構影響されてる気もしますが…今回はこの121倍って結果になりました。

実際に計算速度とは関係ないのですが、C言語はプログラムの実行に1分30秒程、Pythonの方はプログラムを事項するのに1時間ほどかかりました。(なぜかtimeコマンドの結果と違う。。。)
これの原因としては、Pythonはファイルから読み込んだ数値を格納する際に動的に配列サイズを大きくしながら確保しています。そこらへんで、ファイルの読み込みに時間がかかってるのかな〜と…。(先に1億個分の配列を確保しようとしたら、こちらもうまく動かなかった。プログラムのミスの可能性がある。。。)

 

おわりに

前回記事で混乱させてしまった方、本当にごめんなさいm(_ _)m。
前回のプログラムではC言語の方が54倍速いって結果になりましたが、今回は160倍です。これは、Pythonで1年間かかる計算がC言語だと3日で終わる計算です。と、言ってもPythonで単純にforループを回すようなプログラムは遅いってのは用がないです。また、数値計算用のライブラリNumPyとか使ったりするとPythonでももっと早くなるのかもしれません。

 

いやぁ〜、プログラミング楽しい✌(‘ω’✌ )三✌(‘ω’)✌三( ✌’ω’)✌

 

訂正

• Pythonのプログラムでファイルから読み込んだ値が使われていない
• Pythonの時間計測をCと同じtime.clock()を使用したほうがいいのでは？

とのコメントを頂きました。これらについて訂正しました。
いやぁ〜、Pythonでファイルから読み込んだ値を使ってないってのはヒドイ…。なぜ自分で気づけなかったのか…。具体的には42行目から52行目あたりを修正しました。
部分的に修正前と修正後で表示したいところですが、直す箇所が多すぎるのでそのままブログを書き換えて更新しました。

ご指摘ありがとうございますm(_ _)m

訂正箇所について、プログラムを訂正して計測しなおした(ファイルから読み込んだ値を使う)プログラムの方が、計算時間が長くなるという結果となりました。これは、配列から値を取ってくる動作の方が乱数を生成するよりも時間がかかっているという事になります。(1億個も配列にいれてるから？)

ほげぇ〜まじか〜。