Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
binomial modulo (verified: codechef, SANDWICH)
#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));
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment