[bzoj 3738] [Ontak2013]Kapitał

求C(n+m, n)去掉所有末尾的0后对10^k取模的结果. (1≤n,m≤10^15, 1≤k≤9) 组合数对素数的幂取模的一种方法见 [bzoj 2142] 礼物.

题外话: 以前在hzwer的博客看到一个对阶乘质因数分解的方法 (黄学长又说来自一位叫怡红公子的高人), <具体数学> 上也看到了, 和上述组合数取模的方法是一致的.

仍然考虑中国剩余定理. 怎么去掉末尾所有的0呢? 用上述方法求出, 末尾0的数目是. 只需在中消去多于的2或5, 补上不足的2或5:

此处我是分类讨论的......学习了Po姐的代码才发觉可以更优雅.

最后把模的答案和模的答案合并即可.

#include <bits/stdc++.h>
#define For(i, a, b) for (int i = (a), i##_ = (b); i <= i##_; ++i)

using namespace std;
typedef long long ll;
const int P[] = {2, 5}, N = 1953126;
ll f[2][N], M[3];

void init(ll n)
{
	For (i, 0, 1) {
		f[i][0] = 1;
		For (j, 1, n)
			f[i][j] = f[i][j-1] * (j % P[i] ? j : 1) % M[i];
	}
}

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

inline ll inverse(ll x, ll n)
{
	ll y, t;
	exgcd(x, n, y, t); // yx + tn = 1
	return (y + n) % n;
}

ll fpm(ll x, ll n, ll m)
{
	if (n < 0) {
		x = inverse(x, m);
		n = -n;
	}
	ll y = 1;
	for (; n; n >>= 1, (x *= x) %= m)
		if (n & 1)
			(y *= x) %= m;
	return y;
}

ll fac(ll n, ll& e, ll* const &f, const ll& p, const ll& m)
{
	if (n == 0) return 1;
	e += n/p;
	return fpm(f[m], n/m, m) * f[n%m] % m * fac(n/p, e, f, p, m) % m;
}

inline ll C(ll n, ll m, int i, ll& e)
{
	e = 0;
	ll t = 0, a = fac(n, e, f[i], P[i], M[i]),
		b = fac(m, t, f[i], P[i], M[i]) * fac(n-m, t, f[i], P[i], M[i]) % M[i];
	e -= t;
	return a * inverse(b, M[i]) % M[i];
}

int main()
{
	int k;
	ll n, m;
	scanf("%lld%lld%d", &n, &m, &k);
	M[0] = 1<<k;
	M[1] = 1;
	For (i, 1, k) M[1] *= 5;
	M[2] = M[0] * M[1];

	init(M[1]);

	ll a, b, x = C(n+m, n, 0, a), y = C(n+m, n, 1, b);
	if (a > b) {
		(x *= fpm(2, a-b, M[0]) * fpm(5, -b, M[0]) % M[0]) %= M[0];
		(y *= fpm(2, -b, M[1])) %= M[1];
	} else {
		(x *= fpm(5, -a, M[0])) %= M[0];
		(y *= fpm(5, b-a, M[1]) * fpm(2, -a, M[1]) % M[1]) %= M[1];
	}

	ll ans = M[1] * inverse(M[1], M[0]) % M[2] * x % M[2]
		+ M[0] * inverse(M[0], M[1]) % M[2] * y % M[2];
	ans -= ans >= M[2] ? M[2] : 0;

	printf("%0*lld\n", k, ans);
	return 0;
}