2つの自然数が互いに素である確率にπが!?【プログラムで実験して確かめる】python

記事の内容

今回の記事では、数学的にもプログラミング的にもとても興味深いテーマを紹介します。

数学youtuberの鈴木貫太郎さんが紹介していた話です。

www.youtube.com

それは、

2つの自然数が互いに素である確率は、 $\dfrac{6}{\pi ^{2}}$ である、というものです。

もう少し説明しましょう。
自然数は無限にあります。その中から無作為に一つの自然数を選ぶ。さらに、重複を許しつつもう一回自然数を無作為に選びます。

このとき選ばれた2つの自然数が互いに素である確率を考えるのです。確率の問題として考えるということに注意しましょう。

だからこそ、円周率である $\pi$ が登場することが不思議に見えてしまいます。

くわえて、無限は厄介です。無限があるせいでイメージし辛くなります。

具体的に考えてみましょう。

範囲を100までの自然数とします。

この範囲で自然数を一つ選ぶ。さらに、もう一回選ぶ。

5と19が選ばれたとします。このとき、この2つの自然数が互いに素になるかどうかを考えるのです。

5と19なら、互いに素です。
21と63なら、互いに素ではありません。
68と7なら、互いに素です。

このように試行を繰り返せば、そこには確率が現れてきます。

少しはイメージがわいてきたでしょうか？

それでは目次をごらんください。

記事の内容
どこからが出てくるの？
プログラムなら、実験できるのでは！
Pythonによる実装
関連記事

どこから $\pi$ が出てくるの？

それでは、もう少し数学的な話に入っていきましょう。

どこから $\pi$ が出てくるのか？
この謎の答えは、「自然数の平方の逆数の和」のせいです。あの天才オイラーが解決した問題です。
オイラーは導いた答えは次の数式です。

$\dfrac{1}{1^{2}}+\dfrac{1}{2^{2}}+\dfrac{1}{3^{2}}+\dfrac{1}{4^{2}}+\ldots =\dfrac{\pi ^{2}}{6}$

「自然数の平方の逆数の和」を考えるために、三角関数で近似しているのです。三角関数を利用しているのならば、円周率であるπが登場するのも自然です。

詳しい照明は、鈴木貫太郎さんが説明してくれています。どのように三角関数を利用しているのか、ぜひチェックを。

www.youtube.com

そして、2つの自然数が互いに素である確率を求めようとすると、「自然数の平方の逆数の和」が式に現れてくるのです。よって、答えに $\pi$ が出てくることになります。

プログラムなら、実験できるのでは！

ここで一つ思いついたことがあります。
確率の話なのだから、実際に繰り返し試行してやればいい。数式の結果を確かめられるはずです。繰り返すことが得意な奴といったらコンピュータの出番です。コンピュータでも無限は扱えませんが、桁を増やしてあげれば近似することができます。

コンピュータに繰り返し試行させることで、実際に確率を求めてみましょう！

そのために、今回はpythonで実装してみました。

ポイントは互いに素であるかどうかの判定です。これは、互いに1以外の公約数を持たなければいいのです。

2数を与えると最大公約数を返してくれる関数を考えます。このとき最大公約数が1ならば、互いに素であると判定できます。最大公約数を求めるアルゴリズムといえば、ユークリッドの互除法です。

ユークリッドの互除法では、以下の手順で最大公約数を求めます。

x を y で割り、余り r を求める
y を r で割り、その余り r1 を求める
r を r1 で割り、その余り r2 を求める

（余りが0になるまで繰り返す）

余りが 0 になったときの割る数が最大公約数になる

再帰で実装するのがよさそうですね。

Pythonによる実装

import math #円周率をつかうため
import random

def gcd1(a1, b1): #ユークリッドの互除法による最大公約数判定
    if b1 == 0:
        return a1
    else:
        return gcd1(b1, a1% b1)

s = 0
num = 1

while num <= 10000000:
  a = random.randint(1,10000000)
  b = random.randint(1,10000000)
  if gcd1(a,b) == 1: #最大公約数が1のとき互いに素であると判定
    s += 1
  num += 1


print(s / num) #実験値
print(6 / math.pi**2) #理論値