mattyuuの数学ネタ集

世界は数式で溢れている

素数大富豪〜会心の一撃の成功確率〜

はじめに

素数大富豪が数学好きの間で静かなブームになっています。「素数大富豪って何ですか?」という方は、本アドベントカレンダー1日目のにせいさんの記事を読んでください。

さて、今年もアドベントカレンダーの季節がやってきました! この記事は、素数大富豪アドベントカレンダー 1日目の記事です。

素数大富豪アドベントカレンダー! - にせいの日記


私は素数大富豪を1回しかやったことがなく、戦略やゲームの中で出会った素数達について語れるレベルではありませんが、せっかくの機会ですので何か投稿したいと思い、素数大富豪 AdventCalendar 2016に2回分登録させていただきました。本記事はその1つ目の投稿です。2つ目は12月24日です。

では、よろしくお願いします。

問題

突然ですが、素数大富豪には一発大逆転の秘策があり、例えば対戦相手の手札が1枚、自分の手札が10枚で圧倒的に不利な状況でも、うまくヒットすれば逆転勝利が可能です。これは通常の大富豪にはない素数大富豪の特徴です。

その秘策とは"秘"密にするほど大したものではなく、自分からカードを出せる場面で、手持ちの全カードを一列に並べて場に出すことです。場に出したカードから決まる数が素数であれば、手札が全てなくなるため勝利となります。しかし合成数だった場合は、ペナルティとして山札から出したカードと同じ枚数だけ引かないといけないためハイリスク、ハイリターンの策といえるでしょう。

ここでは、自分の手札全てを一気に場に出すこの策を「会心の一撃」と呼ぶことにします(上述の通り会心の一撃に失敗すると大ダメージをくらうため、メガンテと呼んだ方が妥当な気もします。)。

先日のMathPowerで(第一回)素数大富豪大会が開かれましたが、MathPowerの司会も務められたタカタ先生が会心の一撃を狙って場を盛り上げていた記憶があります。今回の記事では、会心の一撃について考察します。タイトル通りですが、問題は下記です。

f:id:mattyuu:20161206014435j:plain

答えの確率をNを使った簡単な数式で表すことができれば大変よいのですが、ただでさえ不規則な素数が考察対象になる上に、トランプという自然界から見たらなんの必然性もないよく分からない意味不明で超人工的な存在(数は1~13の13種類で、数字は人間の指が10本あるためか10進法で、スートというよく分からない記号があるため各数のカードの枚数は4)も考察対象になることから、それは不可能に思えます。今回はコンピュータを用いて頑張ることにします。

シミュレーションによる解法

N枚のトランプで会心の一撃を行うときの成功確率は下記で求めることができます。


\displaystyle
\frac{\mbox{N枚のカードを一列に並べた時に素数となるパターン数}}{\mbox{N枚のカードを一列に並べる全パターン数}}


さて、まずは分母の全パターンを計算してみましょう。
それは手札にN枚のカードが配られていて、それを一列に並べるパターンとなります。
まず52枚のカードからN枚を選出するパターンが{}_{52}{\rm C}_N通り、その各々に対してN枚のカードを一列に並べるパターンはN!通りになるので、

\displaystyle
{}_{52}{\rm C}_{N}*N!=\frac{52\cdot 51 \cdots (52-N+1)}{N!}\cdot N!=52\cdot 51\cdots (52-N+1) \bigl(={}_{52}{\rm P}_N\bigr)

となります。(最初から{}_{52}{\rm P}_Nで計算すればよかったですね。)

{}_{52}{\rm P}_Nのパターンの数のうち、何パターンが素数になるのかをコンピュータに計算させることで確率の分子を得ることができます。そこから手札がN枚の時の確率を求めますることができるのですが、この分母の{}_{52}{\rm P}_N、つまり素数判定しないといけないパターンの数が、Nと共にどのくらいの数になるのかというと、

 N=1:52\mbox{通り}
 N=5:311,875,200\mbox{通り}
 N=10:57,407,703,889,536,000\mbox{通り}
 N=20:306,532,583,223,109,543,994,300,006,400,000\mbox{通り}
 N=30:71,759,895,859,162,404,648,138,410,532,134,505,676,800,000,000\mbox{通り}
 N=40:168,388,112,212,869,181,588,664,081,406,834,062,715,634,990,448,640,000,000,000\mbox{通り}
 N=50:40,329,087,585,471,939,285,830,318,428,201,883,487,644,752,720,441,638,912,000,000,000,000\mbox{通り}
 N=52:80,658,175,170,943,878,571,660,636,856,403,766,975,289,505,440,883,277,824,000,000,000,000\mbox{通り}

となります。N=52では8,065不可思議通りもの素数判定をしなければいけません。そして、その判定する桁数は68桁です(1桁*9(1~9)*4(スート)+2桁*4(10~K)*4(スート))。そんな巨大な数を素数判定するだけでも時間的なコストがかかりそうです。。。

実際N枚のカードで全通りのパターンに対して素数判定を行う方法では、N=8(←うろ覚え)で私のコンピュータと私の我慢の限界に達しました。
1日かかっても終わりません。N=20くらいまでいくと、現在のスーパーコンピュータでも宇宙誕生から今日までの日数をかけても計算できないかもしれません。

そこで、今回は各Nに対してそれぞれ10万回ずつ会心の一撃を繰り出し、その成功確率を算出しました。(実際は10万回のシミュレーションを10セット行っていますので、各Nに対して100万回シミュレーションしています。)

コンピュータを使って問題を解かないといけない時点で悔しいですが、その上正確な値ではなく数値実験の結果に止まることは甚だ遺憾です。しかし、世の中は解ける問題よりも解けない問題の方がはるかに多いのです。あきらめるしかないでしょう。

一つ気がかりがあります。N=52では8,065不可思議通りものパターンがあります。そんなにパターンが多いのにたった10万回の施行でそれなりに正しい確率は出せるのでしょうか? 10万回1セットのシミュレーションの結果の確率が上下すれば、それは真の確率に近いものとは言えないでしょう。

この辺りは専門的な確率統計の知識があれば定量的に評価できるものと思われますが、実際に10万回やれば結果の確率にぶれは無く、確率はある値に近づいていくことがわかりました。関東地方の視聴率は関東地方1600万世帯のうち600世帯だけでしか調査していないとのことなので、同じ感覚で10万回やれば十分なのでしょう!

それにしても回数が全然違いますが、、、この辺りは今後統計検定の資格を狙って勉強していきたいと思います。
80,658,175,170,943,878,571,660,636,856,403,766,975,289,505,440,883,277,824,000,000,000,000
100,000

結果

ではシミュレーション結果をお見せします。各Nに対して、10万回の会心の一撃を繰り出しその成功回数を表にしております。上述の通り10万回1セットだけでは本当に有効なシミュレーションができているかわからないため10セット繰り返しています。一番右の列が各Nに対してトータル100万回の施行での確率です。

f:id:mattyuu:20161205220325p:plain

N会心の一撃が成功した回数はセットによらずおおよそ一定となっており、それなりに有効なシミュレーションになっていそうということはわかっていただけると思います。Nに対して確率がどのように変動するのかがわかるグラフも書いておきます。

f:id:mattyuu:20161206014449j:plain

見ていただければわかるかと思いますが、Nが大きくなるにつれて会心の一撃の成功確率は急激に下がっていきます。
N=1では成功確率46%(13種類の数のうち、素数は2,3,5,7,11,13の6種類なので6/13≒46%)ですが、N=10で成功確率は4%に落ちます。
素数定理によると、整数x以下の素数の個数\pi(x)は、

\displaystyle
\pi(x)\simeq\frac{x}{\log x}

で表せるので、整数x素数である確率p(x)

\displaystyle
p(x)\simeq\frac{\pi(x)}{x}=\frac{1}{\log x}

のように考えることができます。xが大きくなればなるほどp(x)はどんどん小さくなります。つまり整数が大きくなっていくにつれて素数はまばらになっていくのです。Nを大きくするということはより大きな整数で会心の一撃を行うわけで、出した数が素数になる確率はどんどん下がっていくのです。このp(x)と今回得られた確率のグラフと比べるのは乱暴ですが、似ているような気がしませんか?

f:id:mattyuu:20161206014508j:plain

さて、得られた結果をよく見ると、N=52会心の一撃の成功確率が上がっています。これは一体どうしたことでしょうか? 答え(と予想するもの)は注釈に書きますので、みなさん考えてみてください。*1

曲がりなりにも答えの確率を求めることができましたが、不満に思われている方もいるでしょう。素数大富豪に限らず、2の倍数、3の倍数、5の倍数は簡単に判定できます。会心の一撃を繰り出す時に、その判定を行わない人はいないでしょう。

ということで、会心の一撃で出したトランプから決まる数が2、3、5のいずれかの倍数の場合は、ノーカウントととしてシミュレーションしたグラフも載せておきましょう。

f:id:mattyuu:20161206014524j:plain

グラフだけでは分かりにくいですが、N=52でも確率は2.5%程度までしか落ちません。40回に1回は成功する計算になります。

まとめ、感想

冒頭に記載した、素数大富豪アドベントカレンダー5日目のにせいさん記事を見ると、上級者は会心の一撃を狙ってないことがわかり、少し意外でした。確かに今回の結果を見ても会心の一撃の成功率は高くありません。大会となると人の目も気になるので、会心の一撃を出しにくいのかもしれませんね。

ただ、素数大富豪は究極の素数出会い系カードゲームですので、今後街中で会心の一撃に失敗した人を見かけたら、笑ったり、哀れんだりせずに、「全力で未知の素数に出会おうとしている肉食系の人だ」と心の中で褒めてあげましょう。

おまけ

最後に今回コンピュータを使ってどのようにシミュレーションしたのかを、超簡単にお伝えします。

私は普段Javaを使っているのですが、Java素数大富豪で出てくるような巨大な数を扱うのは難しいと思っています。BigIntegerというクラスを使わないといけないですし、BigIntegerに素数判定のために用意されているメソッドも大きい整数に対しては確率的にしか判定できないようです。(当初BigInteger用の素数判定メソッドを作りかけていたのですが、絶対遅いだろうと思って途中でやめました。)

そのため、以前から少し触っていたSageというアプリを使うことにしました。インストールしたのが昔ですので、どのサイトを参照してどうインストールしたかは覚えていません。インストールしたい方は、「Sage 数学」で検索してみてください。

Sageを使えば\zeta(2)1000!も簡単に求めることがきます。

f:id:mattyuu:20161206015221p:plain

is_prime()というメソッドが用意されており、次のような巨大な数でも素数判定を一瞬で行ってくれます。

f:id:mattyuu:20161206015205p:plain

Sageでプログラミングをする時はPythonというプログラミング言語を使うようです。今回初めてPythonでプログラミングしました。今回のシミュレーションのソースは下記の通りです。

import random

# n枚のトランプでkaisuu回会心の一撃を繰り出して、その成功回数を出力する
def sosuDaifugo(n,kaisuu):
	# 成功した回数を格納する変数
	counter = 0

	# kaisuu回繰り返す
	for i in range(0,kaisuu):
		# トランプの定義
		trump = [ "1", "1", "1", "1", "2", "2", "2", "2", "3", "3", "3", "3", "4", "4", "4", "4", "5",
				"5", "5", "5", "6", "6", "6", "6", "7", "7", "7", "7", "8", "8", "8", "8", "9", "9", "9", "9", "10", "10",
				"10", "10", "11", "11", "11", "11", "12", "12", "12", "12", "13", "13", "13", "13"]

		num = ""
		for j in range(0,n):
			# ランダムにトランプを配るイメージ
			card = trump[random.randint(0,len(trump)-1)]
			# 配られたトランプを右に並べていくイメージ。文字列結合であることに注意
			num = num  + card
			# 配られたトランプは山から無くすイメージ
			trump.remove(card)

		# 素数判定
		if is_prime(int(num)):
			# ここのコメントを外すと出会えた素数も出力される
			# print int(num)
			counter = counter + 1

	print counter


# メイン部分 1枚から52枚までのループ
for n in range(1,53):
	print n
	# 100000回繰り返す
	sosuDaifugo(n,100000)

このソースファイルをsosu.py (.pyはPythonソースコードの拡張子)として保存し、Sage上でloadすれば実行されます。

f:id:mattyuu:20161206015118p:plain

素数大富豪のシミュレーションと聞くと大変そうなイメージを持つ方が多いかもしれませんが、ソースは意外に短いでしょう?こんな簡単なソースでシミュレーションできるわけですから、素数大富豪に限らず身の回りで気になることがある方はどんどんシミュレーションしていきましょう!

f:id:mattyuu:20161206015636j:plain

*1:「52枚のカードを全て出すと、各桁の和は256となり、3の倍数にならない。つまり52枚全て使ったときの数も3の倍数にならない。3の倍数が出なくなる分、心持ち素数が出やすくなる」と予想します。