bzoj 4830: [Hnoi2017]抛硬币 [范德蒙德卷积 扩展lucas]

时间：2022-10-15 21:56:17

4830: [Hnoi2017]抛硬币

题意：A投a次硬币，B投b次硬币，a比b正面朝上次数多的方案数，模\(10^k\)。

\(b \le a \le b+10000 \le 10^{15}, k \le 9\)

几乎一下午和一晚上杠这道题...中间各种翻《具体数学》各种卡常

有两种做法，这里只说我认为简单的一种。

题目就是要求

\[\sum_{i=0}^a \sum_{j=0}^b [i>j] \binom{a}{i} \binom{b}{j}
\]

化一化得到

\[\sum_{i=1}^a \sum_{j=0}^{a-i} \binom{a}{i+j} \binom{b}{b-j}
\]

根据范德蒙德卷积，得到

\[\sum_{i=b+1}^{a+b} \binom{a+b}{i}
\]

然而并不会计算...根据题目范围显然应该与a,b的差相关...

参考了这里

把这个组合数拆成两部分计算：

\[\sum_{i=\lceil \frac{a+b}{2} \rceil}^{a+b} \binom{a+b}{i}+ \sum_{i=b+1}^{\lfloor \frac{a+b}{2} \rfloor} \binom{a+b}{i}
\]

\(a+b\)为奇数时前一部分就是\(2^{a+b-1}\)。

注意判断\(a+b\)为偶数的时候，要减去\(\frac{1}{2} \binom{a+b}{(a+b)/2}\)

因为模数特殊需要使用扩展lucas定理

但是\(10^9\)直接用会T...

10只有2和5两个质因子，预处理不能整除2或5的数的阶乘就没问题了。

关于模意义再提两点：

  1. \(a \equiv b \pmod {md} \rightarrow a\equiv b \pmod m\),所以可以放在\(10^9\)下做
  2. 组合数除以2，对于\(2^a\)这个质因子2的逆元不存在，我们只能把指数减1.然后2的指数有可能本来就为0，所以不是所有组合数都能除以2(有的本来就是奇数啊)，可以证明\(\binom{2n}{n}\)为偶数

最后还有一个卡常技巧：

如果\(p\)指数大于9不要继续计算了...节省好多常数

#include <iostream>

#include <cstdio>

#include <cstring>

#include <algorithm>

#include <ctime>

using namespace std;

typedef long long ll;

const int N=2e6+5;

int mo, p2, p5;

ll a, b; int k;

inline ll Pow(ll a, ll b, ll mo) {

	ll ans=1;

	for(; b; b>>=1, a=a*a%mo)

		if(b&1) ans=ans*a%mo;

	return ans;

}

inline ll Inv(ll a, ll n) {return Pow(a, ((n&1) ? n/5*4 : n/2) - 1, n);}

ll fac2[N], fac5[N];

inline ll Fac(ll n, ll p, ll pr) {

	if(n == 0) return 1;

	ll ans = 1;

	if(p == 2) ans = fac2[pr] %pr; else ans = fac5[pr] %pr;

	ans = Pow(ans, n/pr, pr);

	if(p == 2) ans = ans * fac2[n%pr] %pr; else ans = ans * fac5[n%pr] %pr;

	return ans * Fac(n/p, p, pr) %pr;

}

ll cou(ll n, ll p) { ll ans = 0; for(ll i=n; i; i/=p) ans += i/p; return ans; }

ll C(ll n, ll m, bool div2=0) {

	if(n < m) return 0;

	ll ans = 0, a = 0, t;

	{

		ll p = 2, pr = p2;

		ll c = cou(n, p) - cou(m, p) - cou(n-m, p);

		if(div2) t = Pow(p, --c, pr);

		else t = Pow(p, c, pr);

		if(c<9) {

			a = Fac(n, p, pr) * Inv(Fac(m, p, pr), pr) %pr * Inv(Fac(n-m, p, pr), pr) %pr * t %pr;

			a = a * (mo/pr) %mo * Inv(mo/pr, pr) %mo;

			ans = (ans + a) %mo;

		}

	}

	{

		ll p = 5, pr = p5;

		ll c = cou(n, p) - cou(m, p) - cou(n-m, p);

		t = Pow(p, c, pr);

		if(div2) t = t * Inv(2, pr) %pr;

		if(c<9) {

			a = Fac(n, p, pr) * Inv(Fac(m, p, pr), pr) %pr * Inv(Fac(n-m, p, pr), pr) %pr * t %pr;

			a = a * (mo/pr) %mo * Inv(mo/pr, pr) %mo;

			ans = (ans + a) %mo;

		}

	}

	return ans;

}

char c[50];

void print(ll ans) {

	for(int i=0; i<20; i++) c[i] = '0';

	int p=0;

	while(ans) c[++p] = ans%10 + '0', ans /= 10;

	for(int i=k; i>=1; i--) putchar(c[i]); puts("");

}

void init(int k) {

	fac2[0] = 1; p2 = Pow(2, k, 1e9);

	for(int i=1; i<=p2; i++) fac2[i] = fac2[i-1] * (i%2 ? i : 1) %p2;

	fac5[0] = 1; p5 = Pow(5, k, 1e9);

	for(int i=1; i<=p5; i++) fac5[i] = fac5[i-1] * (i%5 ? i : 1) %p5;

}

int main() {

	freopen("in", "r", stdin);

	init(9); mo = 1e9;

	while(scanf("%lld %lld %d", &a, &b, &k) != EOF) {

		ll ans = Pow(2, a+b-1, mo);

		for(ll i=b+1; i<=(a+b)/2; i++) ans = (ans + C(a+b, i)) %mo;

		if(~(a+b)&1) ans = (ans - C(a+b, (a+b)>>1, 1) + mo) %mo;

		int t = 1; for(int i=1; i<=k; i++) t *= 10; ans %= t;

		print(ans);

	}

	//printf("time %lf\n", (double) clock()/CLOCKS_PER_SEC);

}

相关文章