码迷,mamicode.com
首页 > 其他好文 > 详细

[BZOJ 3129] [Sdoi2013] 方程 【容斥+组合数取模+中国剩余定理】

时间:2015-04-06 20:13:58      阅读:383      评论:0      收藏:0      [点我收藏+]

标签:

题目链接:BZOJ - 3129

 

题目分析

使用隔板法的思想,如果没有任何限制条件,那么方案数就是 C(m - 1, n - 1)。

如果有一个限制条件是 xi >= Ai ,那么我们就可以将 m 减去 Ai - 1 ,相当于将这一部分固定分给 xi,就转化为无限制的情况了。

如果有一些限制条件是 xi <= Ai 呢?直接来求就不行了,但是注意到这样的限制不超过 8 个,我们可以使用容斥原理来求。

考虑容斥:考虑哪些限制条件被违反了,也就是说,有哪些限制为 xi <= Ai 却是 xi > Ai,这样就转化为了 xi >= Ai 的限制条件。

那么我们就可以在 2^8 * T(求组合数) 的时间内求出答案了。

怎么求这个组合数呢?直接预处理阶乘的逆元是不可以的,因为模数不都是质数。

我们要将模数拆成一个个 pi^ai 这样的形式,使得它们两两之间互质,就可以分别求出答案,最后再用中国剩余定理组合起来。

中国剩余定理:如果有n个方程 x = xi (mod mi) ,M = m1 * m2 * .. * mn ,那么在 mod M 的意义下,方程组有一个唯一解。

x = sigma(Mi * Inv(Mi) * xi) % M ,其中 Mi = M / mi ,Inv(Mi)是Mi在mod mi意义下的逆元。

那么我们的问题就是,如何求出 C(n, m) % (p^a) 。

这里就需要用到“组合数取模”了,专门用来求解这种问题。

使用类似于快速阶乘的方法,将组合数中分数线上下的阶乘都拆成 e * p^f 的形式,然后 e 直接计算,f 分数线上下相减之后再计算。

怎么将 x! 拆成 e * p^f 呢?

假设我们要 mod 的数是 p^a ,那么我们需要预处理出 [1, p^a - 1] 中除去 p 的倍数的其余数的前缀积(类似阶乘少了 p 的倍数)。

然后我们知道 [1, x] 中包含 p 的数有 x / p 个,我们将这些数中都提取出 1 个 p,那么就获得了 p^(x/p),然后这 x / p 个数就变成了 [1, x/p],就可以递归下去。

其余的部分可以分段来求,分成 [1, p^a - 1], [p^a + 1, p^a + p^a - 1] ..... 这样,每一段的积都是一样的 (mod p^a 意义下),直接快速幂就可以了。

最后还会剩下一段 [1, x % (p^a)] ,也是直接预处理出的值。

这样这道题就做完了(呼~)。

 

另外注意的是,在写代码的时候,我求逆元使用欧拉定理但是确用错了。

欧拉定理:a^phi(b) = 1 (mod b) 条件:gcd(a, b) = 1

注意是 a^phi(b) 而不是 a^(b-1) !当 b 不是质数的时候就跪了!  

 

代码

#include <iostream>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdio>

using namespace std;

typedef long long LL;
typedef double LF;

const int MaxP = 10201 + 15, MaxN1 = 8 + 5;

int T, p, n, n1, n2, m, Top, Ans;
int A[MaxN1];

LL Temp;
LL Fac[10][MaxP], Pr[10], Pi[10], Pa[10], Phi_Pi[10], Mi[10], Inv_Mi[10], Xi[10];

LL Pow(LL a, LL b, LL Mod)
{
	LL ret, f;
	ret = 1; f = a;
	while (b)
	{
		if (b & 1) 
		{
			ret *= f;
			ret %= Mod;
		}
		b >>= 1;
		f *= f;
		f %= Mod;
	}
	return ret;
}

void Prepare()
{
	int x, SqrtX;
	x = p; SqrtX = (int)sqrt((LF)x);
	Top = 0;
	for (int i = 2; i <= SqrtX; ++i)
	{
		if (x % i != 0) continue;
		Pr[++Top] = i;
		Pa[Top] = 0;
		Pi[Top] = 1;
		while (x % i == 0)
		{
			++Pa[Top];
			Pi[Top] *= i;
			x /= i;
		}
		Phi_Pi[Top] = Pi[Top] / Pr[Top] * (Pr[Top] - 1);
	}
	if (x > 1) 
	{
		Pr[++Top] = x;
		Pa[Top] = 1;
		Pi[Top] = x;
		Phi_Pi[Top] = Pi[Top] - 1;
	}
	for (int i = 1; i <= Top; ++i)
	{
		Mi[i] = p / Pi[i];
		Inv_Mi[i] = Pow(Mi[i], Phi_Pi[i] - 1, Pi[i]);
		Fac[i][0] = 1;
		for (int j = 1; j < Pi[i]; ++j) 
		{
			if (j % Pr[i] != 0) Fac[i][j] = Fac[i][j - 1] * j % Pi[i];
			else Fac[i][j] = Fac[i][j - 1];
		}
	}
}

struct ES
{
	LL e, f;
};

ES Calc(int x, int k)
{
	ES ret, tc;
	if (x < Pr[k]) 
	{
		ret.e = Fac[k][x];
		ret.f = 0;
		return ret;
	}
	ret.f = x / Pr[k];
	tc = Calc(x / Pr[k], k);
	ret.f += tc.f;
	ret.e = tc.e * Fac[k][x % Pi[k]] % Pi[k];
	ret.e = ret.e * Pow(Fac[k][Pi[k] - 1], x / Pi[k], Pi[k]) % Pi[k];
	return ret;
}

LL C(int x, int y, int k)
{
	LL ret;
	int pf;
	ES Ex, Ey, Exy;
	Ex = Calc(x, k);
	Ey = Calc(y, k);
	Exy = Calc(x - y, k);
	ret = Ex.e * Pow(Ey.e, Phi_Pi[k] - 1, Pi[k]) % Pi[k] * Pow(Exy.e, Phi_Pi[k] - 1, Pi[k]) % Pi[k];
	pf = Ex.f - Ey.f - Exy.f;
	if (pf >= Pa[k]) ret = 0;
	else ret = ret * Pow(Pr[k], pf, Pi[k]) % Pi[k];
	return ret;			
}

int C(int x, int y)
{
	if (x == y) return 1;
	if (x < y) return 0;
	if (y == 0) return 1;
	LL ret = 0;
	for (int i = 1; i <= Top; ++i) Xi[i] = C(x, y, i);
	for (int i = 1; i <= Top; ++i)
	{
		ret += Xi[i] * Mi[i] % p * Inv_Mi[i] % p;
		ret %= p;
	}
	return (int)ret;
}

void DFS(int x, int Cnt, int Sum)
{
	if (x == n1)
	{
		if (Cnt & 1) Temp -= C(m - Sum - 1, n - 1);
		else Temp += C(m - Sum - 1, n - 1); 
		Temp = (Temp % p + p) % p;
		return; 
	}
	DFS(x + 1, Cnt, Sum);
	DFS(x + 1, Cnt + 1, Sum + A[x + 1]);
}

int Solve()
{
	Temp = 0;
	DFS(1, 0, 0);
	DFS(1, 1, A[1]);	
	return (int)Temp;	
}

int main()
{
	scanf("%d%d", &T, &p);
	Prepare();
	for (int Case = 1; Case <= T; ++Case)
	{
		scanf("%d%d%d%d", &n, &n1, &n2, &m);
		for (int i = 1; i <= n1; ++i) scanf("%d", &A[i]);
		int Num;
		for (int i = 1; i <= n2; ++i)
		{
			scanf("%d", &Num);
			m -= Num - 1;
		}
		if (n1 > 0) Ans = Solve();
		else Ans = C(m - 1, n - 1);
		printf("%d\n", Ans);
	}
	return 0;
}

  

[BZOJ 3129] [Sdoi2013] 方程 【容斥+组合数取模+中国剩余定理】

标签:

原文地址:http://www.cnblogs.com/JoeFan/p/4396432.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!