Python電卓を使ってみよう 第7回 Python電卓でモンテカルロ法とコラッツ予想

今回は円周率を求めるプログラムとコラッツ予想（コラッツ問題）のプログラムです。コラッツ予想については以下のページを参照してください。NHKのEテレの数学番組でも取り上げられたことがあるので聞いたことがあるかもしれません。数学の未解決問題の1つだそうです。
コラッツの問題
https://ja.wikipedia.org/wiki/コラッツの問題
○円周率をどうやって求めるか
　円周率をどうやって求めるか、より効率よく求めるには（計算するには）どうしたら良いのか？というのは古くからある問題です。数学関係の書籍やWebで検索すると実に様々な方法が紹介されています。その中で「そんな方法で求められるの？」というのがモンテカルロ法です。
　モンテカルロ法は地面に描かれた四角い紙の中に描かれた円を用意します。その四角い紙の中に砂つぶを投げて、円の中と外に入った砂つぶの数から円周率を求めます。1粒だけだと誤差が大きすぎますが、砂つぶを大量に投げることでだんだんと円周率に近い値になります。
モンテカルロ法
https://ja.wikipedia.org/wiki/モンテカルロ法
　これをプログラムで行う場合には乱数を使います。乱数でX,Y座標を求め、その座標が円の中にあるかどうかを調べます。座標の原点を中心とする円の方程式は以下のようになります。（一般的にはx*xはxの2乗で示されますが、ここはプログラムということで）
x*x + y*y = r*r
　半径1未満の円であればX,Yの計算は1未満になります。1未満なら円の中、1以上なら円の外という具合に判定します。
　実際のプログラムは以下のようになります。
from casioplot import *
from random import *
clear_screen()
n=10000
count=0
for i in range(n):
x=random()
y=random()
c=(0,0,255)
if (x*x+y*y)<1:
count=count+1
c=(255,64,0)
set_pixel(int(x*190),int(y*190),c)
show_screen()
p=4*count/n
draw_string(200,5,str(p))
show_screen()
　これまでと同様にパソコン側でファイルを作成しておき、電卓に転送します。ファイル名はM1.pyとしてあります。実行すると図のようになります。赤い部分が円の中、青が円の外を示します。第1象限のみなので円の1/4だけが描かれます。数学座標系とコンピューター座標では上下が逆なので図のようになります。
　上記のプログラムは1万回繰り返しますが、任意の数だけ繰り返した方がいいかもしれません。このような場合はinput()を使ってプログラム実行時にユーザーに回数を入力してもらうようにします。
　実際のプログラムは以下のようになります。一ヶ所しか変更していませんので説明は不要でしょう。
from casioplot import *
from random import *
clear_screen()
n=int(input('n='))
count=0
for i in range(n):
x=random()
y=random()
c=(0,0,255)
if (x*x+y*y)<1:
count=count+1
c=(255,64,0)
set_pixel(int(x*190),int(y*190),c)
show_screen()
p=4*count/n
draw_string(200,5,str(p))
show_screen()
　これまでと同様にパソコン側でファイルを作成しておき、電卓に転送します。ファイル名はM2.pyとしてあります。実行すると図のようになります。最初に回数を入力します。回数となる数値を入れたらEXEキーを押します。
○円形で表示する
　円の1/4しか表示されないのは物足りないか、なんか違うなと思った人がいるかもしれません。そこで円形に表示するように改良してみましょう。円形で表示するために円の座標計算の範囲を変更します。先ほどは0〜1.0未満でしたが、今度は-1.0〜1.0未満の範囲にします。この場合、random()ではなくuniform()を使います。()の中に範囲を指定します。
　まず、簡単なプログラムを作って結果を確認します。以下のプログラムは6回ランダムな小数値をシェルに表示します。
from random import *
for i in range(6):
n=uniform(-1,1)
print(n)
　プログラムはパソコンで作成しておきM3.pyという名前のファイルにして保存し電卓に転送します。実行すると以下のようになります。ちなみにパソコンでも動作チェック可能です。
乱数の動作を確認したら座標範囲を変更します。random()をuniform()に変更するだけでなく、もう1ヶ所変更する部分があります。それは点を表示する座標です。電卓の画面サイズは384×192です。正円にするので192×192の座標範囲になります。このため円の中心座標は(96,96)になります。これを計算した座標に加える必要があります。
from casioplot import *
from random import *
clear_screen()
n=10000
count=0
for i in range(n):
x=uniform(-1,1)
y=uniform(-1,1)
c=(0,0,255)
if (x*x+y*y)<1:
count=count+1
c=(255,64,0)
set_pixel(int(x*96)+96,int(y*96)+96,c)
show_screen()
p=4*count/n
draw_string(200,5,str(p))
show_screen()
　プログラムはパソコンで作成しておきM4.pyという名前のファイルにして保存し電卓に転送します。実行すると以下のようになります。
　nの値が固定されていますので、任意の数を入力して計算するように変更します。これはn=10000の部分をinput()にするだけです。最初のプログラムを改良したのと同じなので説明は不要でしょう。
from casioplot import *
from random import *
clear_screen()
n=int(input('n='))
count=0
for i in range(n):
x=uniform(-1,1)
y=uniform(-1,1)
c=(0,0,255)
if (x*x+y*y)<1:
count=count+1
c=(255,64,0)
set_pixel(int(x*96)+96,int(y*96)+96,c)
show_screen()
p=4*count/n
draw_string(200,5,str(p))
show_screen()
プログラムはパソコンで作成しておきM5.pyという名前のファイルにして保存し電卓に転送します。実行すると以下のようになります。
○コラッツ予想
　次にコラッツ予想のプログラムを作りましょう。コラッツ予想の計算ルールは以下のようになります。数値nは正数であることが前提です。
nが偶数なら2で割る
nが奇数なら3をかけて1を足す
　これを繰り返すと最終的に1になるというのがコラッツ予想になります。ルールが簡単なのでプログラムも簡単です。
　実際のプログラムは以下のようになります。
n=int(input('n='))
while n!=1:
if n%2 == 0:
n=n/2
else:
n=n*3+1
print(int(n))
プログラムはパソコンで作成しておきM6.pyという名前のファイルにして保存し電卓に転送します。実行すると以下のようになります。
.