15.7.2 階乗の計算

階乗は二つのアルゴリズムを組み合わせて計算しています。 考え方としては: 階乗の基数部分を計算する；最終段階で2のべき乗の乗法をシフト演算で実行する，というものになります。

nが小さい時には，n!の奇数因子の計算は次のように行えることが分かります。 即ち，nより小さい全ての正の奇数因子の積と， [n/2]!の奇数因子を乗じたものが等しくなります。ここで[x]はxの整数部を意味します。 この関係を繰り返して使用して階乗の計算を行います。下記にその例を示します。

23! = (23.21.19.17.15.13.11.9.7.5.3)(11.9.7.5.3)(5.3)2^{19}

現在の実装では，全ての因数を一本のリストに集め，再帰的ではないループでその積をバンバン計算していきます。

nが大きい時には，素数櫛を通して計算を行います。この時にはPeter Luschnyが提唱したヘルパー関数を 使っています。

                            n
                          -----
               n!          | |   L(p,n)
msf(n) = -------------- =  | |  p
          [n/2]!^2.2^k     p=3

ここでpは素数です。指数kは奇数になるようにしておきます。 というのも，kは[n/2]の2進表現における“1"ビットの数にしておきたいからです。 関数L(p,n)は，pが合成数である時にはゼロとします。任意の素数 pに対しては次のように計算を行います。

          ---
           \    n
L(p,n) =   /  [---] mod 2   <=  log (n) .
          ---  p^i                p
          i>0

このヘルパー関数を使うと，n!=[n/2]!^2*msf(n)*2^kという漸化式を使って，n!の奇数部分の計算ができるようになります。この漸化式は，[n/2^i]上の小さいnに対するアルゴリズムを用いて停止させます。

上記の二つのアルゴリズムはどちらも2進分割法を使ってたくさんの小さい因数の積を求めています。最初に できる限りたくさんの数の積を一つのレジスタに求めておき，マシンワードに収まる因数のリストを生成していきます。 このリストは2分割して，その積を再帰的に求めていきます。

このような分割法を使用することで，繰り返しNx1の積を行うよりは効率的に実行できます。 大きな数の積を作っていくと，Karatsubaアルゴリズムや高次のアルゴリズムを使うことで効率が良くなるからです。 Karatsuba法の閾値より小さい場合は筆算乗算アルゴリズムを使うことになりますが，効率アップを狙って大きなブロック単位で実行するようにできます。