HNOI2017 拋硬幣（組合數學+擴充套件盧卡斯）

題目連結

題目大意

給定$a,b$，考慮所有長度為$a$的01串和長度為$b$的01串，統計所有前者中的1嚴格比後者多的方案數。
答案模$10^k$。
$b\le a\le 10^{15},k\le 9,a-b\le 10,000$

題解

神仙思路……
$a=b$的情況很簡單，只需要考慮有多少組“平局”即可，剩下的勝敗可以一一對應。
平局顯然就是$\sum_{i=1}^a\binom ai\binom bi=\binom{2a}{a}$，因此答案為$\frac{2^{2a}-\binom{2a}{a}}{2}$。
考慮$a>b$，此時我們會發現，任意一種不勝的情況，把所有位異或1必然對應一種勝的情況。於是此時需要考慮的僅僅是有多少局勝的全部異或1也仍然勝利（神仙！）。
令$i$表示$b$中1的個數，$i+j$表示$a$中1的個數，則需要滿足：$j>0,b-i<a-i-j$，化簡可得$0<j<a-b$。
因此，勝轉勝的總方案數如下：
$p=\sum_{i=0}^b\sum_{j=1}^{a-b-1}\binom{a}{i+j}\binom{b}{i}=\sum_{j=1}^{a-b-1}\sum_{i=0}^b\binom{a}{i+j}\binom{b}{b-i}=\sum_{j=1}^{a-b-1}\binom{a+b}{j+b}$
於是就可以在$O(a-b)$的複雜度內計算了。答案就是$\frac{2^{a+b}+p}{2}$。
但是這道題很毒瘤，模數不是質數……我們只能用擴充套件盧卡斯……
更毒瘤的是這道題還卡常，我們算和的時候可以先算出第一個組合數，然後用逆元之類的東西推出後面的，終於AC了……

#include <bits/stdc++.h>
namespace IOStream {
	const int MAXR = 10000000;
	char _READ_[MAXR], _PRINT_[MAXR];
	int _READ_POS_, _PRINT_POS_, _READ_LEN_;
	inline char readc() {
	#ifndef ONLINE_JUDGE
		return getchar();
	#endif
		if (!_READ_POS_) _READ_LEN_ = fread(_READ_, 1, MAXR, stdin);
		char c = _READ_[_READ_POS_++];
		if (_READ_POS_ == MAXR) _READ_POS_ = 0;
		if (_READ_POS_ > _READ_LEN_) return 0;
		return c;
	}
	template<typename T> inline void read(T &x) {
		x = 0; register int flag = 1, c;
		while (((c = readc()) < '0' || c > '9') && c != '-');
		if (c == '-') flag = -1; else x = c - '0';
		while ((c = readc()) >= '0' && c <= '9') x = x * 10 - '0' + c;
		x *= flag;
	}
	template<typename T1, typename ...T2> inline void read(T1 &a, T2&... x) {
		read(a), read(x...);
	}
	inline int reads(char *s) {
		register int len = 0, c;
		while (isspace(c = readc()) || !c);
		s[len++] = c;
		while (!isspace(c = readc()) && c) s[len++] = c;
		s[len] = 0;
		return len;
	}
	inline void ioflush() { fwrite(_PRINT_, 1, _PRINT_POS_, stdout), _PRINT_POS_ = 0; fflush(stdout); }
	inline void printc(char c) {
		_PRINT_[_PRINT_POS_++] = c;
		if (_PRINT_POS_ == MAXR) ioflush();
	}
	inline void prints(char *s) {
		for (int i = 0; s[i]; i++) printc(s[i]);
	}
	template<typename T> inline void print(T x, char c = '\n') {
		if (x < 0) printc('-'), x = -x;
		if (x) {
			static char sta[20];
			register int tp = 0;
			for (; x; x /= 10) sta[tp++] = x % 10 + '0';
			while (tp > 0) printc(sta[--tp]);
		} else printc('0');
		printc(c);
	}
	template<typename T1, typename ...T2> inline void print(T1 x, T2... y) {
		print(x, ' '), print(y...);
	}
}
using namespace IOStream;
using namespace std;
typedef long long ll;
typedef pair<int, int> P;

const int MAXN = 2000005, MAXM = 10005;
ll pf2[MAXM], pf5[MAXN], aa[MAXM], bb[MAXM];
void extgcd(ll a, ll b, ll &x, ll &y) {
	if (b > 0) {
		extgcd(b, a % b, y, x);
		y -= a / b * x;
	} else x = 1, y = 0;
}
ll modpow(ll a, ll b, int c) {
	ll res = 1;
	for (; b; b >>= 1) {
		if (b & 1) res = res * a % c;
		a = a * a % c;
	}
	return res;
}
int get_inv(int x, int p) {
	ll v, y;
	extgcd(x, p, v, y);
	return v;
}
int calc_fact(ll *pf, ll n, int p, int q) {
	ll res = 1, cnt = 0;
	for (; n; n /= q) {
		res = res * pf[n % p] % p;
		if (n >= p) cnt += n / p;
	}
	return res * modpow(pf[p], cnt, p) % p;
}
ll calc_cnt(ll n, int p) {
	ll res = 0;
	for (n /= p; n; n /= p) res += n;
	return res;
}
int calc_binom(ll *pf, ll n, ll m, int p, int pk) {
	if (m > n) return 0;
	int a = calc_fact(pf, n, pk, p), b = calc_fact(pf, m, pk, p), c = calc_fact(pf, n - m, pk, p);
	ll t = calc_cnt(n, p) - calc_cnt(m, p) - calc_cnt(n - m, p);
	return (ll)a * get_inv(b, pk) % pk * get_inv(c, pk) % pk * modpow(p, t, pk) % pk;
}
ll solve(ll *pf, ll a, ll b, int p, int pk) {
	if (a == b + 1) return 0;
	int ta = calc_fact(pf, a + b, pk, p), tb = calc_fact(pf, a - 1, pk, p), tc = calc_fact(pf, b + 1, pk, p);
	ll cnt = calc_cnt(a + b, p) - calc_cnt(a - 1, p) - calc_cnt(b + 1, p);
	ll t = (ll)ta * get_inv(tb, pk) % pk * get_inv(tc, pk) % pk, res = t * modpow(p, cnt, pk) % pk;
	for (int i = 2; i < a - b; i++) {
		ll x = a - i + 1, y = b + i;
		for (; x % p == 0; x /= p) ++cnt;
		for (; y % p == 0; y /= p) --cnt;
		t = t * (x % pk) % pk * get_inv(y % pk, pk) % pk;
		res += t * modpow(p, cnt, pk) % pk;
	}
	return (res % pk + pk) % pk;
}
int main() {
	ll a, b; int k;
	while (~scanf("%lld%lld%d", &a, &b, &k)) {
		int p2 = modpow(2, k + 1, 1000000000), p5 = modpow(5, k, 1000000000), pp = p2 * p5;
		for (int i = pf2[0] = 1; i <= p2; i++) {
			pf2[i] = pf2[i - 1];
			if (i % 2) pf2[i] = pf2[i] * i % p2;
		}
		for (int i = pf5[0] = 1; i <= p5; i++) {
			pf5[i] = pf5[i - 1];
			if (i % 5) pf5[i] = pf5[i] * i % p5;
		}
		int inv1 = get_inv(p5, p2), inv2 = get_inv(p2, p5);
		char prt[10] = "%00lld\n"; prt[2] = '0' + k;
		if (a == b) {
			ll x = calc_binom(pf2, a + b, b, 2, p2);
			ll y = calc_binom(pf5, a + b, b, 5, p5);
			x = inv1 * x % p2 * p5;
			y = inv2 * y % p5 * p2;
			x = ((x + y) % pp + pp) % pp;
			printf(prt, (modpow(2, a + b - 1, pp) - x / 2 + pp) % (pp / 2));
		} else {
			ll x = solve(pf2, a, b, 2, p2), y = solve(pf5, a, b, 5, p5);
			x = inv1 * x % p2 * p5, y = inv2 * y % p5 * p2;
			x = ((x + y) % pp + pp) % pp;
			printf(prt, (modpow(2, a + b - 1, pp) + x / 2) % (pp / 2));
		}
	}
	return 0;
}