Python電卓でモンテカルロ法とコラッツ予想

コラッツの問題
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という名前のファイルにして保存し電卓に転送します。実行すると以下のようになります。

• 

• 

• 

　計算結果が1行ごとに表示されてしまうので途中経過が確認しにくい状態です。改行せずに,で区切って表示するように変更します。pythonではprint()の()内にend=','とすることで末尾に,をつけることができます。
　改良したプログラムは以下のようになります。

n=int(input('n='))
while n!=1:
  if n%2 == 0:
    n=n/2
  else:
    n=n*3+1
  print(int(n),end=',')
print()

プログラムはパソコンで作成しておきM7.pyという名前のファイルにして保存し電卓に転送します。実行すると以下のようになります。

• 

• 

• 

　いろいろな数を入力して1に収束するかやってみてください。と書きたいところですが、fx-CG50では入力した数値によってはメモリ不足なのかエラーで停止してしまうことがあります。このプログラムはパソコンでも動くので電卓でエラーだったらパソコンで動かしてみてください。パソコンであれば問題なく動きます。

• 

• 

• 

著者　仲村三郎
日本のお家芸である小型化を体現したのが古き良き時代のポケットコンピューター（略してポケコンと呼ばれていました）。関数電卓はポケットコンピューターとは違いますが、今ではPythonでプログラミングできるまでに進化。進化した関数電卓を使って楽しもうではありませんか。