Decomposition of a factorial to a * 2^k mod 2^64

factorial.cc

#include <bits/stdc++.h>
using int64 = unsigned long long;

typedef std::vector<int64> poly;

const int M = 1024;

int64 comb[64][64], ok[M];
poly operator*(const poly& A, const poly& B) {
  poly C(64);
  for (int i = 0; i < 64; ++i) {
    if (A[i]) {
      for (int j = 0; i + j < 64; ++j) {
        C[i + j] += A[i] * B[j];
      }
    }
  }
  return C;
}
poly shift(const poly& A, int64 B) {  // f(x) -> f(x + B)
  std::vector<int64> P(64);
  for (int i = 0; i < 64; ++i) {
    P[i] = i == 0 ? 1 : P[i - 1] * B;
  }
  poly C(64);
  for (int i = 0; i < 64; ++i) {
    if (A[i]) {
      for (int j = 0; j <= i; ++j) {
        C[j] += A[i] * comb[i][j] * P[i - j];
      }
    }
  }
  return C;
}

poly power(int64 x) {
  poly C(64);
  if (x == 0) {
    C[0] = 1;
    return C;
  } else {
    poly A = power(x / 2);
    poly B = A * shift(A, x / 2);
    if (x % 2 == 1) {
      C[0] = 2 * x - 1;
      C[1] = 2;
      return B * C;
    } else {
      return B;
    }
  }
}

int64 power(int64 A, int64 B) {
  int64 res = 1;
  while (B) {
    if (B & 1) {
      res = 1LL * res * A;
    }
    A = 1LL * A * A;
    B >>= 1;
  }
  return res;
}

std::unordered_map<int64, int64> h;

// \Prod_{i=1}^{x}{2i - 1} mod 2^64 => 
// \Prod_{i=1}^{x}{2x - 1} mod x^64
int64 dp(int64 x) {
  if (h.count(x)) {
    return h[x];
  }
  return h[x] = power(x)[0];
}

// return {f, s} as n! = 2^f * s mod 2^64
std::pair<int64, int64> fun(int64 x) {
  if (x == 0) {
    return {0, 1};
  }
  int64 odd = dp((x + 1) / 2);
  std::pair<int64, int64> wtf = fun(x / 2);
  wtf.first += x / 2;
  wtf.second *= odd;
  return wtf;
}

void work() {
  int64 n;
  scanf("%lld", &n);
  printf("%lld\n", dp(n));
  std::pair<int64, int64> a = fun(n + n);
  std::pair<int64, int64> b = fun(n);

  if (a.first - b.first - (n - 1) >= 64) {
    puts("0");
  } else {
    int64 res = a.second;
    res <<= a.first - b.first - (n - 1);
    // printf("%lld %lld\n", res, b.second);
    res *= power(b.second, (1ULL << 63) - 1);
    printf("%lld\n", res);
    for (int i = 63; i >= 0; --i) {
      putchar('0' + (res >> i & 1));
    }
    puts("");
  }
}

int main() {
  int T = 1;
  for (int i = 0; i < 64; ++i) {
    for (int j = 0; j <= i; ++j) {
      comb[i][j] = j == 0 ? 1 : comb[i - 1][j - 1] + comb[i - 1][j];
    }
  }

  scanf("%d", &T);
  while (T--) {
    work();
  }
  return 0;
}