LOJ#2023. 「AHOI / HNOI2017」抛硬币（OGF+ExLucas+Crt）

https://loj.ac/problem/2023

题解：

考虑生成函数做：
(sum_{y>0} (1+x)^a*(1+x^{-1})^b [x^y])
(=sum_{y>b} (1+x)^{a+b} [x^y])

注意到(a-b)非常小，也就是说(b)靠近((a+b)/2)

由组合数对称性，我们知道((sum_{i=0}^{n/2} inom{n}{i})=frac{2^n+[n~is~even]*inom{n}{n/2}}{2})

那么就只用多算(5000)项组合数了，注意(2^9)和(5^9)很小，直接套用(crt+ExLucas)即可。

(ExLucas)需要卡常。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

void exgcd(ll a, ll b, ll &x, ll &y) {
	if(b == 0) { x = a, y = 0; return;}
	exgcd(b, a % b, y, x); y -= x * (a / b);
}

ll inv(ll a, ll p) {
	ll x, y;
	exgcd(a, -p, x, y);
	x = (x % p + p) % p;
	return x;
}

ll a, b, k;

const int N = 2e6 + 5;

template<int p, int mo>
struct nod {
	ll ksm(ll x, ll y) {
		ll s = 1;
		for(; y; y /= 2, x = x * x % mo)
			if(y & 1) s = s * x % mo;
		return s;
	}
	
	struct P {
		ll x, y;
		P(ll _x = 0, ll _y = 0)	{
			x = _x, y = _y;
		}
	};

	P mtp(P a, P b) {
		return P(a.x * b.x % mo, a.y + b.y);
	}
	
	ll fac[N];
	
	void build() {
		fac[0] = 1;
		fo(i, 1, mo) fac[i] = fac[i - 1] * ((i % p == 0) ? 1 : i) % mo;
	}
	
	P jc(ll n) {
		if(n < p) return P(fac[n], 0);
		ll x = fac[n % mo];
		if(fac[mo] == mo - 1 && (n / mo) % 2 == 1) x = x * (mo - 1) % mo;
		P w = jc(n / p);
		w.y += n / p; w.x = w.x * x % mo;
		return w;
	}
	
	P fan(P a) { return P(inv(a.x, mo), -a.y);}
	
	ll C(ll n, ll m) {
		if(n < m) return 0;
		P a = mtp(mtp(jc(n), fan(jc(n - m))), fan(jc(m)));
		if(a.y >= k) return 0;
		return a.x * ksm(p, a.y) % mo;
	}
	
	ll C2(ll n, ll m) {
		if(n < m) return 0;
		P a = mtp(mtp(jc(n), fan(jc(n - m))), fan(jc(m)));
		a.y --;
		if(a.y >= k) return 0;
		return a.x * ksm(p, a.y) % mo;
	}
	
	ll solve() {
		ll s = ksm(2, a + b - 1);
		if((a + b) % 2 == 0) {
			if(p == 2) s -= C2(a + b, (a + b) / 2); else
				s -= C(a + b, (a + b) / 2) * inv(2, mo) % mo;
		}
		for(ll x = b + 1; x <= (a + b) / 2; x ++) s = (s + C(a + b, x)) % mo;
		
		s = (s % mo + mo) % mo;
		return s;
	}
};

nod<2, 512> t0;
nod<5, 1953125> t1;

ll m1, c1, m2, c2;

ll gcd(ll x, ll y) {
	return !y ? x : gcd(y, x % y);
}

void merge() {
	ll t = gcd(m1, m2);
	m2 /= t;
	if((c2 - c1) % t != 0) return;
	ll a = inv(m1 / t % m2, m2) * ((c2 - c1) / t % m2 + m2) % m2;
	c1 = a * m1 + c1;
	m1 = m1 * m2;
}

int main() {
	freopen("coin.in", "r", stdin);
	freopen("coin.out", "w", stdout);
	t0.build(); t1.build();
	while(scanf("%lld %lld %lld", &a, &b, &k) != EOF) {
		m1 = 1; fo(i, 1, k) m1 *= 2;
		m2 = 1; fo(i, 1, k) m2 *= 5;
		c1 = t0.solve() % m1;
		c2 = t1.solve() % m2;
		merge();
		int cw = 0;
		for(ll x = c1; x; x /= 10) cw ++;
		cw = max(cw, 1);
		fo(i, 1, k - cw) pp("0");
		pp("%lld", c1);
		hh;
	}
}