Raspberry Pi スパコン (6) 円周率


2023年 02月 16日

乱数を用いた円周率の計算

複数のコンピュータ(ノード)で乱数を用いた円周率の計算を行い、すべての結果の平均をとれば、より正確な円周率が求まるはず。

次の図のように、1辺が1の正方形を用意し、左下隅を原点(0,0)とする。
この原点を中心とする半径1の円弧をこの正方形の上に描くと、正方形は、1/4円の扇型(図の赤の部分)と、残り(図の青の部分)に分けられる。
この正方形の上にランダムに点をとると、その点が赤い部分に含まれる割合はπ/4であるから、その割合を4倍すれば円周率になる。

プログラムは簡単なので、説明は省略する。

# -*- coding: utf-8 -*-
import random
import matplotlib.pyplot as plt

def near_pi(n):
	inside = 0

	for i in range(n):
    	x = random.random()
    	y = random.random()

    	if x*x+y*y < 1.0:
       	inside += 1
       	plt.plot(x, y,"r.",markersize=2)
    	else:
       	plt.plot(x, y,"b.",markersize=2)
    
	return 4.0 * inside / n

if __name__ == '__main__':
	print( "円周率: {:.10f}".format(near_pi(10000)) )

	# 結果を表示
	plt.gca().set_aspect('equal', adjustable='box')

	plt.grid(True)
	plt.show()

実行すると、しばらくして円周率を表示し、上の図を表示する。

点の数は、1万個にしてある。

$ python3 near_pi.py
円周率: 3.1352000000

MPIを使って、RP0で2つのノードを動かしたら、こうなった。

$ mpiexec -H RP0,RP0 python3 near_pi.py
円周率: 3.1456000000
円周率: 3.1496000000

マスターのRP0以外のワーカーマシンを指定すると、エラーになってしまう。

端末を開いてコマンドを実行したマシンRP0以外でpyplotを動かそうとすると、応答がなくなってしまうようだ。

$ mpiexec -H RP0,RP1 python3 near_pi.py
Traceback (most recent call last):
  File "near_pi.py", line 3, in <module>
	import matplotlib.pyplot as plt
ModuleNotFoundError: No module named 'matplotlib'

ということで、pyplotを使わず、円周率の近似値だけ出力するプログラムに変更してみた。
なお、1万回のループを1億回にし、ホスト名とrankも表示するようにした。

# -*- coding: utf-8 -*-

from mpi4py import MPI
import socket
import random

name = socket.gethostname()
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

def near_pi(n):
	inside = 0

	for i in range(n):
    	x = random.random()
    	y = random.random()
    	if x*x+y*y < 1.0:
       	inside += 1
    
	return 4.0 * inside / n

if __name__ == '__main__':
	print( "{} rank {:2d}   円周率:{:.10f}".format(name,rank,near_pi(100000000)) )
$ time mpiexec -H RP0,RP1,RP2,RP3,RP4,RP5,RP6,RP7 python3 near_pi2.py
RP5 rank  5   円周率:3.1417414800
RP6 rank  6   円周率:3.1415374000
RP4 rank  4   円周率:3.1414466800
RP7 rank  7   円周率:3.1414066400
RP2 rank  2   円周率:3.1415328800
RP3 rank  3   円周率:3.1414802800
RP0 rank  0   円周率:3.1414867600
RP1 rank  1   円周率:3.1417481600

real    1m37.616s
user    1m36.757s
sys    0m0.191s

複数のノードで別々に円周率の近似値を求めることができているようだ。
1億回で100秒なので、1回の処理に1マイクロ秒の時間が掛かっている。
円周率の制度だが、3.141までは値が一致している。
これらの近似値を、マスター(RP0)に集めて平均を出せば、さらに真の円周率に近づくのではないだろうか。でも、各ノードでprintしたものが、マスターの端末に表示されているだけで、まだマスターは各ワーカーの結果を受け取っていない。

データの授受は次回に説明しよう。