binomial modulo (verified: codechef, SANDWICH)

binomial_mod.cc

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <cstdio>
#include <algorithm>
#include <functional>

using namespace std;

#define fst first
#define snd second
#define all(c) ((c).begin()), ((c).end())
#define TEST(s) if (!(s)) { cout << __LINE__ << " " << #s << endl; exit(-1); }

typedef long long ll;


ll add(ll a, ll b, ll M) { // a + b (mod M)
  return (a += b) >= M ? a - M : a; 
}
ll sub(ll a, ll b, ll M) { // a - b (mod M)
  return (a -= b) <  0   ? a + M : a; 
}
ll mul(ll a, ll b, ll M) { // a * b (mod M)
  ll r = a*b - (ll)((long double)(a)*b/M+.5)*M;
  return r < 0 ? r + M: r;
}
ll div(ll a, ll b, ll M) { // solve b x == a (mod M)
  ll u = 1, x = 0, s = b, t = M; 
  while (s) { // extgcd for b x + M s = t
    ll q = t / s;
    swap(x -= u * q, u);
    swap(t -= s * q, s);
  }
  if (a % t) return -1; // infeasible
  return mul(x < 0 ? x + M : x, a / t, M); // b (xa/t) == a (mod M)
}
ll pow(ll a, ll b, ll M) { // a^b mod M
  ll x = 1;
  for (; b > 0; b >>= 1) {
    if (b & 1) x = mul(x, a, M);
    a = mul(a, a, M);
  }
  return x;
}

// (n, k) mod M. complexity: O(\sum_j Q_j) 
// where Q_j = p_j^{e_j} is a prime power divisor.
ll binomial(ll n, ll k, ll M) {
  auto factorial_order = [&](ll n, ll p) {
    ll e = 0, q = p;
    while (q <= n) {
      e += n / q;
      q *= p;
    }
    return e;
  };
  // (n, k) mod p^a by Granville's formula
  auto binomial_prime_power = [&](ll n, ll k, ll p, ll a) {
    ll r = n - k;
    ll b = factorial_order(n, p)
         - factorial_order(k, p)
         - factorial_order(r, p);
    if (b >= a) return 0ll;
    ll c = a - b, pa = pow(p, a), pc = pow(p, c);
    ll acc = 1;
    vector<ll> fact = {1}; 
    for (ll x = 1; x <= pc; ++x) {
      if (x % p) acc = mul(acc, x, pc);
      fact.push_back(acc);
    }
    ll num = 1, den = 1;
    bool is_negative = false;
    for (int d = 0; n; ++d) {
      if (acc > 0 && d >= c) {
        is_negative ^= (n & 1);
        is_negative ^= (r & 1);
        is_negative ^= (k & 1);
      }
      num = mul(num, fact[n % pc], pc);
      den = mul(den, fact[r % pc], pc);
      den = mul(den, fact[k % pc], pc);
      n /= p;
      r /= p;
      k /= p;
    }
    ll val = div(num, den, pc);
    if ((p != 2 || c < 3) && is_negative) val = pc - val;
    return mul(pow(p, b, pa), val, pa);
  };
  // factorize, solve, and chinese remainder
  ll x = 0;
  vector<ll> a, m;
  for (ll p = 2, N = M; N > 1; ++p) {
    ll e = 0, pe = 1;
    while (N % p == 0) {
      N /= p; pe *= p; ++e;
    }
    if (e > 0) x = add(x, 
                   mul(binomial_prime_power(n, k, p, e), 
                   mul(M/pe, div(1, M/pe, pe), M), M), M);
  }
  return x;
}

// compare
ll binom(ll n, ll k, ll mod) {
  vector<vector<ll>> dp(n+1, vector<ll>(k+1));
  function<ll(ll,ll)> rec = [&](ll n, ll k) {
    if (k == 0 || n == k) return 1ll;
    if (dp[n][k] > 0) return dp[n][k];
    return dp[n][k] = add(rec(n-1,k-1),rec(n-1,k),mod);
  };
  return rec(n, k);
}

int main() {
  int ncase; scanf("%d", &ncase);
  for (int icase = 0; icase < ncase; ++icase) {
    srand(icase);
    ll n, k, m;
    scanf("%lld %lld %lld", &n, &k, &m);
    if (n % k == 0)  printf("%lld %lld\n", n/k, 1ll);
    else if (k == 1) printf("%lld %lld\n", n/k+1, 0ll);
    else printf("%lld %lld\n", n/k+1, binomial(k-(n%k)+(n/k), n/k, m));
  }
}