数列・級数・組み合わせ論
数列とは、$a_0, a_1, a_2, \dots , a_n$のように数が並んだもの。
あるいは無限に続く数列の和である無限級数も考える。
組合せ論では漸化式や確率・統計で使われるnPr, nCrの話もする。
import
import sympy as sym
oo = sym.oo # 無限大
# 使用する変数の定義(小文字1文字は全てシンボルとする)
(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) = sym.symbols('a b c d e f g h i j k l m n o p q r s t u v w x y z')
import scipy as sci
数列
sympyでの数列はsym.sequence(数式, (添字シンボル, 開始, 終了))で生成できる。
問: 数列 $a_n=n^2 (0 \le n \le 10)$を生成せよ。
A = sym.sequence(n**2, (n, 0, 10))
print([int(N) for N in A])
print(A[3])
[0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
9
数列が0～10で規定されているからAのリストのサイズは11。
問: 数列 $b_n=\frac{1}{n}$を生成せよ。
B = sym.sequence(1/n)
print(B[10**40])
1/10000000000000000000000000000000000000000
範囲を指定しないと無限数列になるがイテラブルではない。
要素は配列のようにインデックスで取得できる。
2つの数列の和
C = A+B
C[5]
$$\frac{126}{5}$$
AとBは定義域が異なっていたが普通に足せた。
もちろんC[11]はエラーだ。
級数
数列から級数に
sequenceオブジェクトのformulaプロパティを参照すると一般項に戻せる。
C.formula
$$n^2+\frac{1}{n}$$
和はsym.summation(一般項の式, (シンボル, 開始, 終了))で計算できる。
sym.summation(C.formula,(n,1,k))
$$\frac{1}{6} \left(2 k^{3} + 3 k^{2} + k + 6 \operatorname{harmonic}{\left (k \right )}\right)$$
harmonicは調和級数である。
R = sym.summation(r**i,(i,1,n)) # 初項r, 公比rの等比数列
R
$$\begin{cases} n & \text{for}\: r = 1 \\\frac{r - r^{n + 1}}{- r + 1} & \text{otherwise} \end{cases}$$
otherwiseの部分とかも自動出力してくれる。
あなたが神か！
Σ記号のまま関数的に扱ってみる
Sx = sym.Sum(n**x, (n, 1, k))
Sx
$$\sum_{n=1}^{k} n^{x}$$
sym.factor(Sx.subs([(x, 10)]).doit())
$$\frac{k}{66} \left(k + 1\right) \left(2 k + 1\right) \left(k^{2} + k - 1\right) \left(3 k^{6} + 9 k^{5} + 2 k^{4} - 11 k^{3} + 3 k^{2} + 10 k - 5\right)$$
代入、遅延評価、因数分解を行っている。
無限級数
$$\sum_{n=1}^{\infty} \frac{n}{e^n}$$
sym.summation(n/sym.exp(n), (n, 1, oo))
$$\frac{1}{e \left(- \frac{1}{e} + 1\right)^{2}}$$
交代級数
$$\sum_{n=1}^{\infty} \frac{(-1)^n}{n}$$
sym.summation((-1)**n/n, (n, 1, oo))
$$-\log(2)$$
ゼータ関数
$$\sum_{n=1}^{\infty} \frac{1}{n^s}$$
で定義される無限級数をゼータ関数という。
$s=1$のとき調和級数となる。
sym.Sum(1/n**s, (n, 1, oo))とは本質的には異なる挙動を示します。
sym.zeta(-1)
$$-\frac{1}{12}$$
$1+2+3+\dots=-1/12$が出ました。
組み合わせ論
階乗
$$n!=\prod_1^n n$$
sym.factorial(20)
$2432902008176640000$
n=10000ぐらいならすぐに返してくれる。
ガンマ関数
sym.gamma(21.0000000001)
$2.43290200891149 \cdot 10^{18}$
階乗の概念を実数まで拡張したもの。
$n!=gamma(n+1)$
nPr
n個のものからr個を選んで並べる順列。
sym.functions.combinatorial.numbers.nP(10, 4)
$5040$
nCr
n個のものからr個を選ぶ組み合わせ。
sym.functions.combinatorial.numbers.nC(10, 4)
$210$
一般的なプログラミングシーンでのnPr、nCr
数値計算などではscipyのを使うとよいのでは。
整数か浮動小数の近似でもいいのかで選ぶ。
print(scipy.special.perm(600, 28, exact=True))
print(scipy.special.perm(600, 28))
print(scipy.misc.comb(100000, 10, 1))
print(scipy.misc.comb(100000, 10, 0))
323855237805857758873023874171467769493658465635822219454064291636838400000000
3.23855237806e+77
27544920827561470257469913571105409574990000
2.75449208276e+43
分割数
sym.functions.combinatorial.numbers.nT(7, 3) # 分割数
$4$
7を3つに分割するのは5+1+1, 4+2+1, 3+3+1, 3+2+2の4通りがあるからnT(7, 3)=4
##カタラン数
$$C_n=\frac{(2n)!}{(n+1)!n!}$$
漸化式が$C_0=C_1=1, C_{n+1}=C_0C_n+C_1C_{n-1}+\dots+C_nC_0$となるもの。
[sym.catalan(x) for x in range(15)] # カタラン数
[1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012,742900,2674440]
フィボナッチ数
$$F_{n+2}=F_{n+1}+F_n,\ F_0=0,\ F_1=1$$
の漸化式で定義される数。
[sym.fibonacci(x) for x in range(15)] # フィボナッチ数
[0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
リュカ数
$$F_{n+2}=F_{n+1}+F_n,\ F_0=2,\ F_1=1$$
の漸化式で定義される数。
[sym.lucas(x) for x in range(15)] # リュカ数
[2,1,3,4,7,11,18,29,47,76,123,199,322,521,843]
黄金比
$\phi=\frac{1+\sqrt{5}}{2}$は隣り合うフィボナッチ数の比の極限を示す。
sym.N(sym.GoldenRatio, 50)
$1.6180339887498948482045868343656381177203091798058$