(2n-1)!! (奇数階乗) を mod 2^64 で求める

いろいろ忘れないようにメモ✍しておくとよさそうなので書いておきます。

やること

$n \geq 1$ に対して $2n-1$ 以下の奇数をすべて乗算した奇数階乗 (奇数に対する二重階乗)

\begin{align}
(2n-1)!! = 1 \cdot 3 \cdot 5 \cdot \cdots (2n-3)(2n-1)
\end{align}

を $\mod 2^{64}$ で高速に求めます。具体的には $w=64$ として $O(w^2 \log n)$ 時間です。

ここで $w^2$ の部分は多項式の乗算と平行移動 ( $P(x)$ の係数から $P(x+k)$ の係数を求める) にかかります。

やりかた

ソースコードです。Petrozavodsk Camp の問題で出たんですがジャッジがどこか分からない……。いちおう $n \leq 10^6$ の範囲で愚直なのと比較して verify してあります。

uint64_t oddfact(uint64_t n) {
  if (n == 0) return 1;

  static uint64_t comb64[64][64];
  if (comb64[0][0] == 0) {
    for (int i = 0; i < 64; ++i) {
      comb64[i][0] = comb64[i][i] = 1;
      for (int j = 1; j < i; ++j) {
        comb64[i][j] = comb64[i - 1][j] + comb64[i - 1][j - 1];
      }
    }
  }

  using poly = array<uint64_t, 65>;
  poly p = {1};
  uint64_t k = 0;

  for (int b = 63 - __builtin_clzll(n); b >= 0; --b) {
    if (n & (1ULL << b)) {
      for (int i = 63; i >= 0; --i) {
        p[i + 1] += 2 * p[i];
        p[i] *= 2 * k + 1;
      }
      ++k;
    }
    if (b == 0) break;

    poly q = {}, r = {};
    for (int i = 0; i < 64; ++i) {
      uint64_t pk = 1;
      for (int j = i; j >= 0; --j) {
        q[j] += p[i] * comb64[i][j] * pk;
        pk *= k;
      }
    }
    for (int i = 0; i < 64; ++i) {
      for (int j = 0; i + j < 64; ++j) {
        r[i + j] += p[i] * q[j];
      }
    }
    p = r;
    k *= 2;
  }

  return p[0];
}

やってること

多項式 $P_k(x)$ を以下のように定めます。要は $2x+1$ から始まる $k$ 個の奇数の積です。

\begin{align}
P_k(x) = (2x+1)(2x+3)\cdots(2x+2k-1)
\end{align}

すると求めたい $(2n-1)!!$ は $P_n(0)$ になるので、これを計算すればよいです。ここで、すべての $x$ の出現に $2$ が掛かっていることから、 $\mod 2^{64}$ で考える場合 $x^{64}$ 以上の高次の係数はすべて $0$ になります。したがって $0$ 次から $63$ 次までの係数を持っておけば $P_k(x)$ を評価するには十分です。