階乗と二重階乗
階乗 とは次式で定義される演算です。
\[n!=n\,(n-1)\,(n-2) \cdots 2 \cdot 1, \quad 0!=1\]
たとえば、$5$ の階乗は
\[5!=5\cdot 4\cdot 3\cdot 2\cdot 1=120\]
のように計算します。階乗計算は単純な再帰アルゴリズムで実装できます：
\[f(n)=nf(n-1)\]
たとえば $4!$ を計算するときには、$f=4$ からスタートして、$f$ と自然数 $3,\ 2,\ 1$ の積を順次 $f$ 自身に代入していきます。
\[\begin{align*}&f=4\\[6pt]&f=f\times 3\\[6pt]&f=f\times 2\\[6pt]&f=f\times 1\end{align*}\]
Python で階乗関数を実装してみましょう。
# PYTHON_FACTORIAL
# In[1]
# 階乗関数
def fact(x):
if x == 0:
return 1
else:
return x * fact(x - 1)
print(fact(30))265252859812191058636308480000000
math モジュールには階乗を計算する factorial() が用意されています。
# In[2]
import math
# 30の階乗を計算
x = math.factorial(30)
print(x)265252859812191058636308480000000
二重階乗 は
\[\begin{align*}n!!&=n\,(n-2)\,(n-4) \cdots 3 \cdot 1\quad (n=\mathrm{odd\:number})\\[6pt]n!!&=n\,(n-2)\,(n-4) \cdots 4 \cdot 2\quad (n=\mathrm{even\:number})\\[6pt]0!!&=(-1)!!=1\\[6pt]\end{align*}\]
によって定義されます。たとえば、$5$ の二重階乗を計算すると
\[5!!=5\cdot 3\cdot 1=15\]
となります。
math.factorial()
math.factorial(x) は x の 階乗 を返します。
引数 x に非整数や負数を渡すと ValueError を返します。
# MATH_FACTORIAL
# In[1]
# mathモジュールをインポート
import math
# 10の階乗を計算
a = math.factorial(10)
print(a)3628800
scipy.special.factorial()
scipy.special.factorial() は整数 n を受け取って n の階乗 を返します。
デフォルト (exact = False) では処理速度を優先して近似計算を実行しますが、exact に True を渡すと正確な値を返します。
# SCIPY_FACTORIAL
# In[1]
from scipy import special
n = 50
# n!を近似計算
fact_n_1 = special.factorial(n)
# n!の正確な値を計算
fact_n_2 = special.factorial(n, True)
print(fact_n_1)
print(fact_n_2)30414093201713378043612608166064768844377641568960512000000000000
mpmath.factorial()
mpmath.factorial() あるいは mpmath.fac() を使うと、任意精度で階乗計算を実行できます。
# MPMATH_FACTORIAL
# In[1]
from mpmath import *
# 精度を10桁に設定
mp.dps = 10
# 10の階乗
x = factorial(10)
print(x)3628800.0
引数に巨大数を受け取った場合は、スターリングの近似公式を使って階乗の近似値を計算します。
# In[2]
# スターリングの公式を使って巨大数の階乗を計算
x = factorial(10**15)
print(x)1.178796412e+14565705518096756
mpmath.fac2()
mpmath.fac2(x) は $x$ の二重階乗 $x!!$ を返します。
# PYTHON_FACTORIAL_2
# In[1]
from mpmath import *
# 精度を10桁に設定
mp.dps = 10
# 15!!を計算
x = fac2(15)
print(x)2027025.0
sympy.factorial()
sympy.factorial(x) は x の階乗を返します。
# SYMPY_FACTORIAL
# In[1]
import sympy
# 10!を計算
f = sympy.factorial(10)
print(f)362880
SymPy は記号計算に対応するパッケージなので、引数に文字式を渡すと factorial オブジェクトが生成されます。
# In[2]
# 記号nを定義
sympy.var('n')
# (n+2)!を計算
s = sympy.factorial(n + 2)
print(s)factorial(n + 2)
Jupyter notebook の display() を使うと、LaTeX で表示されます (環境によっては sympy.init_printing が必要)。
# In[3]
sympy.init_printing()
display(s)$(n+2)!$ を $n!$ で割ると $(n+2)(n+1)$ になることを確認してみましょう。
# In[4]
# n!を計算
t = sympy.factorial(n)
# (n+2)!をn!で割って数式を簡略化する
expr = sympy.simplify(s/t)
print(expr)(n + 1)*(n + 2)
sympy.factorial2()
sympy.factorial2(x) は x の 二重階乗 を返します。
# SYMPY_FACTORIAL_2
# In[1]
import sympy
# 10の二重階乗を計算
f = sympy.factorial2(10)
print(f)3840