Pythonの階乗の応用

階乗の活用について興味深いテーマをPythonで確認していきます。

ウイルソンの定理

ウイルソンの定理は、階乗の計算を使って素数を求めようとするものです。一見関係のなさそうな両者を結びつけるすごい定理です。

ウイルソンの定理

$(p-1)!\equiv p-1(mod p)$であればpは素数であり、逆も成り立つ。

1つの例で

1つの素数で確認する

素数であるp=5で試してみます。

  1. p=5
  2. factorial=1
  3. for i in range(2,p):
  4. factorial*=i
  5. print(factorial,factorial%p,p-1) 

24 4 4

2.　階乗の値をfactorialで計算することとし、どんどん掛けていくので初期値として1を代入します。
3.4. 1からp-1まで間の階乗を求めます。range(2,p)とすることで、4!=1（初期値）×2×3×4(rangeなのでp=5の1つ前まで)により計算できます。
5.　階乗は24で、5で割った余りは4となりこれはp-1と一致し、素数であることと辻褄があいます。

合成数についても確認する

合成数は、2つの計算結果が異なることを確認するため、6から8までについても確認します。計算したい値をリストp_lに代入し、1つずつ上記と同じ計算をします。

  1. p_l=[6,7,8]
  2. for p in p_l:
  3. factorial=1
  4. for i in range(1,p):
  5. factorial*=i
  6. print(p,factorial,factorial%p,p-1) 

6 120 0 5
7 720 6 6
8 5040 0 7

6と8は合成数で、階乗をpで割った値は0でp-1とは異なります。一方、素数である7の場合、両者は一致することがわかります。とりあえずうまくいったようですが、複数の値について調べるときに、毎回、階乗を計算するのは効率が良くありません。

ウイルソンの定理を使って素数の一覧を作成する

計算の対象を増やし、2から1000までの整数について素数判定をします。上記のプログラムの計算対象を変更し、階乗の計算を工夫します。

  1. prime_list=[]
  2. factorial=1
  3. for i in range(2,1000):
  4. factorial*=(i-1)
  5. if factorial%i==i-1:
  6. prime_list.append(i)
  7. print(prime_list)

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]

階乗はp!=(p-1)!×pで計算することができるので、factorialを前回の(p-1)までの計算結果にpを掛けるようにし、手間を省きます。

計算結果を確認する

前節のプログラムで正しく計算できているか確認します。sympy関数で1000までの素数のリストを作成し、比較します。

  1. import sympy
  2. prime_l2=list(sympy.primerange(2,1000))
  3. prime_list==prime_l2

True

1000までの範囲では、正しく計算できていることが確認できました。

スターリングの近似式

スターリングの近似式を使って階乗の計算をする方法があります。

スターリングの近似式

$n!\simeq\sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n$

単純にスターリングの近似式を計算

スターリングの近似式による階乗の計算
  1. import math
  2. n=5
  3. stir=math.sqrt(2*math.pi*n)*(n/math.e)**n
  4. fact=math.factorial(n)
  5. print('スターリング：',stir,'階乗：',fact)
  6. n=10
  7. s1=math.sqrt(2*math.pi*n)*(n/math.e)**n
  8. stir=math.sqrt(2*math.pi*n)*(n/math.e)**n
  9. fact=math.factorial(n)
  10. print('スターリング：',stir,'階乗：',fact)

スターリング： 118.0191679575901 階乗： 120
スターリング： 3598695.6187410373 階乗： 3628800

階乗の計算なのに、円周率$\pi$が出てくるところが不思議ですが、5!は120に対し約118、10!は3628800に対し約3598695.6と少し小さな結果になってしまいます。

スターリングの近似式改良版

スターリングの近似式の精度を上げるために工夫された改良版があります。上記の式に少し数字を盛っています。

スターリングの近似式 改良版

$n!\simeq\sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n\left(1+\dfrac{1}{12n}\right)$

  1. n=5
  2. stir=math.sqrt(2*math.pi*n)*(n/math.e)**n*(1+1/(12*n))
  3. fact=math.factorial(n)
  4. print('スターリング：',stir,'階乗：',fact)
  5. n=10
  6. stir=math.sqrt(2*math.pi*n)*(n/math.e)**n*(1+1/(12*n))
  7. fact=math.factorial(n)
  8. print('スターリング：',stir,'階乗：',fact)

スターリング： 119.98615409021659 階乗： 120
スターリング： 3628684.7488972126 階乗： 3628800

かなり近い結果になりました。

スターリングの近似式（改良版）と実際の値の差異

改良版スターリングの近似式の精度を確認します。

  1. fact=2
  2. for n in range(3,20):
  3. stir=math.sqrt(2*math.pi*n)*(n/math.e)**n*(1+1/(12*n))
  4. fact*=n
  5. print(f'{int(round(stir,0)):>20d} {int(round(fact,0)):>20d} {fact-stir:>20.3f} {(fact-stir)/fact:.5%}')

                   6                    6                0.002 0.02789%
                  24                   24                0.004 0.01714%
                 120                  120                0.014 0.01154%
                 720                  720                0.060 0.00828%
                5040                 5040                0.314 0.00623%
               40318                40320                1.955 0.00485%
              362866               362880               14.082 0.00388%
             3628685              3628800              115.251 0.00318%
            39915743             39916800             1056.578 0.00265%
           478990872            479001600            10728.204 0.00224%
          6226901267           6227020800           119533.300 0.00192%
         87176841037          87178291200          1450163.326 0.00166%
       1307655337323        1307674368000         19030677.202 0.00146%
      20922521261136       20922789888000        268626864.246 0.00128%
     355683369492896      355687428096000       4058603104.250 0.00114%
    6402308351315582     6402373705728000      65354412418.000 0.00102%
  121643983024716528   121645100408832000    1117384115472.000 0.00092%

nの値が大きくなるにしたがい、絶対的な差異は大きくなりますが、比率は20!で0.00092%とかなり小さくなり、大体の値を計算するだけであれば問題なく使えそうです。