二項展開(binomial expansion)の展開後の各項の係数(二項係数(binomial coefficient))は
階乗(factorial)を利用して以下の用に書けます。
latexの数式
${}_n C _r \equiv \frac{n!}{r!(n - r)!}$
;; Common Lisp
(defun recursive-factorial (n)
(check-type n (integer 0 *))
(if (zerop n)
1
(* n (recursive-factorial (1- n)))))
(defun factorial (n)
(check-type n (integer 0 *))
(if (zerop n)
1
(loop
:named loop
:for i from 1 to n
:for result = i then (* i result)
:finally (return-from loop result))))
(defvar *memoized-factorial* (make-hash-table))
(setf (gethash 0 *memoized-factorial*) 1)
(setf (gethash 1 *memoized-factorial*) 1)
(defun memoized-factorial (n)
(check-type n (integer 0 *))
(let ((x (gethash n *memoized-factorial*)))
(if x
x
(let ((y (* n (memoized-factorial (1- n)))))
(setf (gethash n *memoized-factorial*) y)
y))))
(defun binomial-coefficient (n r)
(check-type n (integer 0 *))
(check-type r (integer 0 *))
(assert (<= r n))
;; n! / r!(n - r)!
(/ (memoized-factorial n)
(* (memoized-factorial r)
(memoized-factorial (- n r)))))
この式はn個の中からr個を取る組み合わせ(combination)の数を求める式でもあります。
要素の順番も考慮する順列(permutation)の数の場合は、以下のようになります。
latexの数式
$ {}_n P _r \equiv \frac{n!}{(n - r)!}$
(defun permutation (n r)
(check-type n (integer 0 *))
(check-type r (integer 0 *))
(assert (<= r n))
;; n! / (n - r)!
(/ (memoized-factorial n)
(memoized-factorial (- n r))))
多くのプログラミング言語では順列を生成するための機能が用意されているようです。
RubyのArrayクラス(array.cで定義)とC++のstd::next_permutationを使ってみます。
> p [1, 2, 3].permutation(2).to_a
[[1, 2], [1, 3], [2, 1], [2, 3], [3, 1], [3, 2]]
## nPr = 3P2 = 3! / (3-2)! = 3! = 6
> p [1, 2, 3].permutation(2).to_a.size
6
> p [1, 2, 3].combination(2).to_a
[[1, 2], [1, 3], [2, 3]]
## nCr = 3C2 = 3! / 2!(3-2)! = 3! / 2 = 3
> p [1, 2, 3].combination(2).to_a.size
3
#include <algorithm>
#include <cstdio>
void print_array(int arr[3]){
printf("%d, %d, %d\n", arr[0], arr[1], arr[2]);
}
int main(void){
int arr[3] = {1, 2, 3};
std::sort(arr, arr+3);
int count = 0;
do {
print_array(arr);
count++;
} while(std::next_permutation(arr, arr+3));
printf("count = %d\n", count);
return 0;
}
// 1, 2, 3
// 1, 3, 2
// 2, 1, 3
// 2, 3, 1
// 3, 1, 2
// 3, 2, 1
// count = 6
これらのアルゴリズムをCommon Lispで書いてみます。
なお、Common Lispには順列を操作するためのライブラリとして cl-permutation があります。
プログラム中で順列を生成したりしたい場合はわざわざ自作せずありがたく利用させて頂きましょう。
;; from Ruby (array.c)
(defun permute0 (n r p index used vals result)
(dotimes (i n)
(when (zerop (bit used i))
;; 順列のindex番目の要素として元の配列のi番の要素を選択
(setf (svref p index) i)
(if (< index (1- r))
(progn
;; 選択済みの要素に対応するフラグをONにする
(setf (bit used i) 1)
(permute0 n r p (1+ index) used vals result)
;; 選択済みフラグをOFFにする
(setf (bit used i) 0))
(progn
;; 添字の配列から値の配列を作成して結果配列に追加
(vector-push
(loop :for j :across p :collect (elt vals j))
result))))))
(defun permute (seq r)
(let* ((n (length seq))
(result (make-array (permutation n r) :fill-pointer 0 :adjustable t))
(used (make-array n :element-type 'bit :initial-element 0))
(p (make-array r :element-type '(integer 0 *))))
(permute0 n r p 0 used seq result)
result))
> (permute '(1 2 3) 2)
#((1 2) (1 3) (2 1) (2 3) (3 1) (3 2))
;; from C++ (std::next_permutation)
(defun next-permutation (arr)
(let ((len (length arr)))
(unless (or (= len 0) (= len 1))
(loop
:for pos :from (1- len) :downto 1
:when (< (aref arr (1- pos)) (aref arr pos))
:do (progn
(rotatef (aref arr (1- pos))
(aref arr
(position-if
#'(lambda (x) (< (aref arr (1- pos)) x))
arr
:from-end t)))
(setf (subseq arr pos)
(nreverse
(make-array (- len pos)
:displaced-to arr
:displaced-index-offset pos)))
(return-from next-permutation t))))))
> (defmacro do-while (test &body body)
`(loop
:do (progn ,@body)
:while ,test))
> (let ((tmp (vector 0 1 2)))
(do-while (next-permutation tmp)
(print tmp)))
#(0 1 2)
#(0 2 1)
#(1 0 2)
#(1 2 0)
#(2 0 1)
#(2 1 0)
他にも色々なアルゴリズムが存在するようです。
その中のひとつとして、階乗進数 を利用したプログラムを書いてみます。
(defun factoradic-permutation (n width)
(let ((fact (make-array width :initial-element 1))
(mantissa (make-array width :initial-element 0)))
;; 階乗を計算して配列に設定
(loop
:for i :from 1 :below width
:for acc := i :then (* acc i)
:do (setf (aref fact i) acc))
;; 階乗進数の仮数を計算して配列に設定
(loop
:for i :from (1- width) :downto 1
:for acc := n :then (mod acc (aref fact (1+ i)))
:do (setf (aref mantissa i)
(floor acc (aref fact i))))
;; 階乗進数を利用して順列を作成 (配列mantissaを使いまわす)
(let ((tmp (loop :for i :from 0 :below width :collect i)))
(loop
:for i :from (1- width) :downto 0
:for x := (aref mantissa i)
:do (setf (aref mantissa i) (nth x tmp)
tmp (delete (nth x tmp) tmp))))
(nreverse mantissa)))
> (dotimes (n 6)
> (print (factoradic-permutation n 3)))
#(0 1 2)
#(0 2 1)
#(1 0 2)
#(1 2 0)
#(2 0 1)
#(2 1 0)
![$V[km/s] = (\frac{398600[km^3/s^2]}{(6378[km] + H[km])})^\frac{1}{2}$](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi1XHiu_7Ap4H7w0r6ZtOW0hlQOBEx25L38yIz2fxPEGMQrvoR-GnEc60vCyX2FRIjf0YFNJOil9x6CpLLkriyKuyxwMqpVbk6POjtL3rTYA51E_ViJ85fnsiqbmpbpYhwSo_z7faqJefsO/s1600/2012-12-16-3_0a7f05a3d858d0ae8763750a02d80a422ca8ce7c.png)
![$T[s] = \frac{2\pi(6378[km] + H[km])}{V[km/s]}$](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhA6LEG6XuayctNeoCeXoz5s2L4c6Mo0Kv5jx9hVWV1T7nwvBP0F7PP-d3EFZcfrcCIrxsLatLEYD0ZaGDCcIMA5go2S3AewnQkujcxhusLm6PzhZf8dSHfkErFkxlmqawlT8c6utcBrQLz/s1600/2012-12-16-3_bea5ab527fb2b63f46da993d5c154613a0860330.png)
