Int 階乗・組合せアルゴリズム
概要
IntCombinatorics クラスは、階乗 $n!$、二重階乗 $n!!$、二項係数 $\binom{n}{k}$ を
静的メソッドとして提供する。
いずれも多倍長整数 Int を返し、任意の大きさの入力に対して正確な結果を得られる。
階乗は入力 $n$ の大きさに応じて最適なアルゴリズムが自動選択される。 小さな $n$ ではテーブル参照で即座に返し、 中程度の $n$ では単純ループ、 大きな $n$ では篩法とトーナメント乗算を組み合わせた高速アルゴリズムを使う。
API の使い方については Int API リファレンス — 組合せ関数 を参照。
階乗アルゴリズム選択
calx は $n$ の値に応じて以下のように階乗アルゴリズムを切り替える。 閾値は実測ベンチマークに基づいて決定されている。
| $n$ の範囲 | アルゴリズム | 計算量 |
|---|---|---|
| $n \leq 100$ | テーブル参照 | $O(1)$ |
| $101 \leq n \leq 228$ | 単純ループ | $O(n \times M(n))$ |
| $n > 228$ | 篩法 + トーナメント乗算 | $O(n \log n \times M(n))$ |
ここで $M(n)$ は $n!$ の桁数に対応するサイズの乗算コストを表す。 Stirling の近似 $\log(n!) \approx n \log n$ より、$n!$ のワード数はおおむね $n$ に比例して増大する。
テーブル参照
$n = 0$ から $n = 100$ までの $n!$ の値をコンパイル時に事前計算したテーブルを持つ。 この範囲の階乗は $O(1)$ で即座に返される。
$100!$ の大きさは:
$$100! \approx 9.3 \times 10^{157}$$
10 進で約 158 桁であり、64 ビットワードでは約 9 ワードに収まる。 テーブル全体のメモリ消費は極めて小さく、キャッシュに十分載るサイズである。
$n \leq 20$ であれば $n!$ は 64 ビット整数に収まるが、
calx は統一的に Int を返すため、$n \leq 100$ の全範囲でテーブル参照を使う。
単純ループ
$101 \leq n \leq 228$ の範囲では、素朴な逐次乗算で $n!$ を計算する。
$$n! = 1 \times 2 \times 3 \times \cdots \times n$$
各ステップで「大きな多倍長整数 $\times$ 小さな整数 (1 ワード)」の乗算を行う。 この乗算は多倍長同士の乗算より格段に軽く、 $n$ ワードの整数に 1 ワードの値を掛ける操作は $O(n)$ で済む。
ループ全体の計算量は $O(n \times M(n))$ であり、 $M(n)$ はここでは「大きな整数 $\times$ 小さな整数」のコストなので実質 $O(n)$ に近い。 ただし結果が成長するにつれて各乗算のコストも線形に増大する。
$n \leq 228$ では、後述する篩法のオーバーヘッド (素数篩の構築、べき乗計算、トーナメントの再帰) が 単純ループの素朴さに勝てないため、この範囲ではループが最速となる。 閾値 228 はベンチマークによって決定された。
篩法 + トーナメント乗算
$n > 228$ では、素因数分解に基づく篩法 (sieve method) と トーナメント乗算 (binary splitting) を組み合わせた高速アルゴリズムを使う。
素因数分解による表現
$n!$ を素因数分解すると:
$$n! = \prod_{p \leq n,\ p \text{ prime}} p^{v_p(n!)}$$
各素数 $p$ に対する指数 $v_p(n!)$ は Legendre の公式で求まる:
$$v_p(n!) = \sum_{k=1}^{\infty} \left\lfloor \frac{n}{p^k} \right\rfloor$$
この和は $p^k > n$ で打ち切られるため、有限ステップで計算できる。
アルゴリズムの手順
- エラトステネスの篩で $n$ 以下の全素数を列挙する。
- 各素数 $p$ に対して Legendre の公式で指数 $v_p(n!)$ を計算する。
- 各素数べき $p^{v_p(n!)}$ を計算する。
- 得られた素数べきの配列をトーナメント乗算で結合する。
トーナメント乗算 (binary splitting)
単純なループ $((a \times b) \times c) \times d$ では、 積のサイズが片側だけ急速に成長し、乗算が不均衡になる。 トーナメント乗算は配列を半分に分割し、再帰的にペアを掛け合わせる:
$$(a \times b) \times (c \times d)$$
この方法ではオペランドのサイズが常に近い値に保たれるため、 Karatsuba や Toom-Cook などの高速乗算アルゴリズムの効率が最大化される。 不均衡な乗算 (巨大な数 $\times$ 小さな数) は 高速乗算の分割統治が十分に効かないため、 サイズの均衡を保つことが重要である。
篩法が優れる理由
単純ループは $1 \times 2 \times 3 \times \cdots \times n$ を逐次的に計算するが、 篩法は以下の点で優位に立つ:
- 冗長な乗算の回避: 合成数 (例: $6 = 2 \times 3$) を個別に掛ける代わりに、 素因数の指数としてまとめて処理する。
- サイズの均衡: トーナメント乗算によりオペランドのサイズが均等に保たれ、 高速乗算アルゴリズムが最大効率で動作する。
計算量: $O(n \log n \times M(n))$
$M(n)$ が超線形 (Karatsuba 以上) である場合、 トーナメント乗算によるサイズ均衡化の効果が顕著になり、 単純ループの $O(n \times M(n))$ を上回る。
二重階乗
二重階乗 $n!!$ は、$n$ から 2 ずつ減らしながら掛け合わせた積である。
偶数 $n$ の場合:
$$n!! = n \times (n-2) \times (n-4) \times \cdots \times 4 \times 2$$
奇数 $n$ の場合:
$$n!! = n \times (n-2) \times (n-4) \times \cdots \times 3 \times 1$$
たとえば $10!! = 10 \times 8 \times 6 \times 4 \times 2 = 3840$、 $9!! = 9 \times 7 \times 5 \times 3 \times 1 = 945$ である。
二重階乗は単純ループで実装している。 $n!!$ の積は $n!$ より大幅に小さい (因子の数が約半分) ため、 篩法やトーナメント乗算のオーバーヘッドが正当化される状況は稀であり、 単純ループで十分な性能が得られる。
二項係数
二項係数 $\binom{n}{k}$ は以下で定義される:
$$\binom{n}{k} = \frac{n!}{k!\,(n-k)!}$$
しかし、この定義通りに 3 つの階乗を計算して除算すると、 中間値が最終結果よりはるかに大きくなり非効率である。 たとえば $\binom{1000}{500}$ は約 299 桁だが、$1000!$ は約 2568 桁に達する。
乗算公式
代わりに、次の乗算公式を使う:
$$\binom{n}{k} = \frac{n \times (n-1) \times \cdots \times (n-k+1)}{1 \times 2 \times \cdots \times k}$$
これは $k$ 個の分子と $k$ 個の分母の積であり、階乗を経由しない。
対称性の利用
$\binom{n}{k} = \binom{n}{n-k}$ であるため、 $k > n/2$ の場合は $k$ を $n - k$ に置き換えて計算する。 これにより乗算回数が最大で半分に削減される。 たとえば $\binom{100}{95}$ は $\binom{100}{5}$ として計算され、 95 回ではなく 5 回の乗除で済む。
逐次乗除算
計算は以下のように逐次的に行う:
$$r_0 = 1, \quad r_{i+1} = r_i \times \frac{n - i}{i + 1} \quad (i = 0, 1, \ldots, k-1)$$
各ステップの $r_i \times (n - i) / (i + 1)$ は必ず整数になることが保証されている。 これは $r_i = \binom{n}{i}$ が常に整数であるという組合せ論的性質による。 したがって、浮動小数点演算や有理数演算を使わず、 整数の乗算と正確な除算のみで計算が完結する。
設計判断
閾値 228 の根拠
篩法 + トーナメント乗算は、素数篩の構築、各素数のべき乗計算、 トーナメントの再帰呼び出しといったオーバーヘッドを伴う。 $n$ が小さい場合は、これらのオーバーヘッドが 単純ループの逐次乗算のコストを上回る。
$n = 228$ がベンチマーク上の損益分岐点であり、 これを超えると篩法の利点 (冗長乗算の排除とサイズ均衡) が オーバーヘッドを上回る。
トーナメント乗算の重要性
$10000!$ を篩法で素因数分解すると、約 1229 個の素数べき ($10000$ 以下の素数の個数) が得られる。 これらの素数べきのサイズは大きく異なる (例: $2^{9995}$ は巨大だが $9973^1$ は小さい)。 逐次乗算では片側だけが急速に成長して不均衡になるが、 トーナメント乗算はペアリングによりオペランドサイズの均衡を保ち、 高速乗算アルゴリズムの効率を最大化する。