\(n\)番目の素数を返すprime(n)

from sympy import prime
a = prime(100)
print(f"100番目の素数は{a}です｡")

## 100番目の素数は541です｡

一行目でsympyモジュールから関数primeをインポートしています｡

素数判定isprime

与えられた数が素数ならTrueを合成数ならFalseを返す関数です｡

from sympy import isprime
isprime(17)

\(x\)以下の素数の個数primepi

from sympy import primepi
print(primepi(100))

## 25

指定の範囲以内の素数をリストするprimerange(a, b)

\(a\)以上で\(b\)より小さいすべての素数のリストを作成します｡ range(a, b)と同じように使えますがprimerange(a, b)のままではリストとして出力できません｡ リストが欲しい場合にはlist()で囲んでください｡

from sympy import primerange
print(list(primerange(10, 20)))

## [11, 13, 17, 19]

from sympy import primerange
for prime in primerange(20, 40):
  print(f"{prime}は素数です｡")

## 23は素数です｡
## 29は素数です｡
## 31は素数です｡
## 37は素数です｡

sieveを用いるとパソコンのメモリーを使いますが, 素数のリストを保存し計算が早くなります｡

from sympy import primerange, sieve
print(list(sieve.primerange(50, 100)))

## [53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

素因数分解factorint(n)

与えられた自然数\(n\)を素因数分解して辞書型で返します｡ Pythonの辞書型についてはまだ話していませんが, ここではfactorint(n)の結果の読み方を説明します｡

from sympy import factorint
a = factorint(1200)
print(a)

## {2: 4, 3: 1, 5: 2}

\(1200\)を素因数分解すると\(2^4\cdot 3 \cdot 5^2\)です｡ これと上の結果{2: 4, 3: 1, 5: 2}を比べると, この辞書型の読み方が見えてくるのではないでしょうか? まずカンマ,で区切られた3つのブロックがあります｡ 最初に2:4, 次に3: 1, そして5: 2です｡ コロン:の左側の数字が\(1200\)の素因数で右側がその指数を表しています｡

他にもたくさん

上で紹介した以外にもSymPyモジュールにはたくさんの便利な関数が入っているので, 気になる人はぜひSymPyの数論に関するドキュメントを読んでみてくださいね｡

\(x\)以下の素数の個数\(\pi(x)\)を\(x\)を変えて表示してみる

｢素数が無限個あることについての考察｣で, \(x\)以下の素数の個数\(\pi(x)\)が\(x\)の値によってどのように変化するかを見ました｡ そのときに使ったデータをSympyモジュールを用いて再現してみましょう｡

下の解答例を見る前に挑戦してみてくださいね｡

解答例

それぞれの\(x\)の値で\(\pi(x)\)を計算してもいいですが, 同じようなことの繰り返しを避けるために, forループを使ってみましょう｡

from sympy import primepi  # sympyモジュールから関数primepiをインポートする
for i in range(1, 7):
  x = 10 ** i              # xは10の1乗から6までを動く
  p = primepi(x)           # 変数pにprime(x)の結果(x以下の素数の個数)を代入
  print(f"π({x})={p}, {x}以下の素数の割合は{p/x}")

## π(10)=4, 10以下の素数の割合は0.4
## π(100)=25, 100以下の素数の割合は0.25
## π(1000)=168, 1000以下の素数の割合は0.168
## π(10000)=1229, 10000以下の素数の割合は0.1229
## π(100000)=9592, 100000以下の素数の割合は0.09592
## π(1000000)=78498, 1000000以下の素数の割合は0.078498

\(4x+1\)型と\(4x+3\)型の素数の個数の比較

\(3\)や\(7\)のように\(4x+3\)の形に書ける素数(\(4x+3\)型)が無限個存在することを｢素数が無限個あることについての考察｣で証明しました｡ 証明方法は違いますが\(4x+1\)型の素数も無限個あることがわかっています｡

では, どちらのグループのほうがより多くの素数を含んでいるのでしょうか? 例えば, \(4x+1\)型の方が\(2\)倍多く素数を含んでいるなど素数分布に偏りはあるか, あるいはどちらのグループにもだいたい同じくらいの素数が入っているのか知りたいですね｡ Pythonを使ってデータを集めてみましょう｡

解答例

\(n\)以下の\(4x+1\)型の素数の個数を示す関数pi_type1を作ってみましょう｡

from sympy import isprime
def pi_type1 (n):
  count = 0
  x = 0
  while 4 * x + 1 <= n:
    if isprime(4 * x + 1):
      count += 1
    x += 1
  return count

\(4x+1\)が与えられた数\(n\)以下になるような\(x\)に対して\(4x+1\)が素数ならcountを1増やすというアルゴリズムです｡

同様に、\(4x+3\)型の素数の個数を示す関数pi_type3を作ってみます｡

from sympy import isprime
def pi_type3 (n):
  count = 0
  x = 0
  while 4 * x + 3 <= n:
    if isprime(4 * x + 3):
      count += 1
    x += 1
  return count

では, これらの関数を使って2つのグループの割合(pi_type1(n)/pi_type3(n))をみてみましょう｡

for i in range(1,7):
  n = 10 ** i
  a = pi_type1(n)
  b = pi_type3(n)
  print(f"{n}以下に4x+1型は{a}個, 4x+3型は{b}個あり, 比は{a/b}です｡")

## 10以下に4x+1型は1個, 4x+3型は2個あり, 比は0.5です｡
## 100以下に4x+1型は11個, 4x+3型は13個あり, 比は0.8461538461538461です｡
## 1000以下に4x+1型は80個, 4x+3型は87個あり, 比は0.9195402298850575です｡
## 10000以下に4x+1型は609個, 4x+3型は619個あり, 比は0.9838449111470113です｡
## 100000以下に4x+1型は4783個, 4x+3型は4808個あり, 比は0.9948003327787022です｡
## 1000000以下に4x+1型は39175個, 4x+3型は39322個あり, 比は0.9962616347083058です｡

この結果を表にしてみました｡

このデータから, 比がほぼ1であることが見て取れます｡ つまり, \(x\)以下の\(4x+1\)型と\(4x+3\)型の素数の個数はだいたい同じくらいありそうだと予測できます｡ 実は, この予測は正しいのです｡ これは算術級数の素数定理で証明されています｡

この定理を証明するためにはまだ数学的準備が足りません｡ この例では, Pythonでのプログラミングを通じて, 素数に関する数多くのデータを簡単に得ることができ, そのデータを使って素数の分布を予想できることを紹介したかったのです｡

科学者が様々な実験でデータを集めるように, 数学者もプログラミングを使ってデータを集め, そこから予想を立て, 実際に厳密な数学的証明を考えることもあります｡

かの有名な数学者カール・フリードリヒ・ガウス(1777ー1855)は一日15分ずつ空いた時間を利用して1000個ずつの自然数にそれぞれいくつの素数が含まれているかを調べ, そのデータに基づいて素数定理を予想したそうです (参考:ウィキペディア:カール・フリードリヒ・ガウス)｡ ガウスの時代にはコンピューターはありませんでしたが, 彼の驚異の計算能力で素数のデータを集めていたんですね｡ 幸い, 私達にはコンピューターが計算を担ってくれるのでデータ集めはすごく手軽になりました｡ このとても便利な道具を使わない手はないですね｡