当前位置：首页 > news >正文

【BJWC2008】王之财宝Gate Of Babylon——超详解

news 2024/10/24 2:45:08

题目描述：P4640 [BJWC2008] 王之财宝

《基尔伽美修》是人类历史上第一部英雄史诗，两河流域最杰出的文学作品之一。作品讲述了基尔伽美修一生的传奇故事。在动画Fate/stay night中，基尔伽美修与亚瑟王等传说中的英雄人物一起出现在了现实世界，展开了一场惊天地、泣鬼神的战斗。在记载于12块泥板的史诗中，基尔伽美修与同伴安吉杜一起降伏了森林的守护者——神兽洪芭芭，成为地上最强的王者，同时将世间所有财宝收归手中。王之财宝(Gate of Babylon)成为Fate中金皮卡（基尔伽美修的外号…）炫耀的资本……

一天金皮卡突发奇想：如果从自己无尽的财宝里面，随便抽不超过M件宝具出来砸死敌人的话。一共有多少种搭配方法呢？假设金皮卡一共有N种不同类型的宝具，大部分类型的宝具都有无限多，但其中T种超级神器的数量是有限的。设第 $i$ 种超级神器的数量不超过 $Bi$ 件。若相同类型的宝具数量相同，则认为是相同的搭配方案。

金皮卡知道方案数会很大，从小数学成绩就好的他挑选了一个质数P，请你帮他计算一下方案数模P后的余数。注意：一件也不选也是一种方案。

输入格式：

第一行包含4个整数，分别为 $N,T,M,P$ 。之后 $T$ 行，每行一个整数，代表 $Bi$ ;

数据范围： $N,M\leq 10^{9},P\leq 10^{5},0\leq Bi\leq 10^{9},M> 0,0\leq T\leq N,T\leq 15.$

输出格式：

仅包含一个整数，即方案数模P后的余数。

分析思路（超详细）：

总体思想是用隔板法+容斥原理求解，具体步骤如下：

1）计算没有任何限制时的所有可能组合数

在不考核任何限制的情况下，N种宝具中每种宝具都可以选择0~M件，总的方案数相当于求解方程： $a+b+c+...+n\leq M$ 的非负整数解的个数，(其中a,b,c...分别代表A、B、C...对应宝具的数量）。可以转化为求解 $a+b+c+...+n+i=M$ (其中虚拟变量 $i$ ，表示多余的“空位”)，这样可以用利用隔板法计算总方案数为： $C^{_{m+n}^{n}}$ 。下面用一个具体示例分析：

比如求共有5个宝具，其中4种超级神器，数量分别为1、2、3、4件，最多拿出10件宝具出来的方案数，即： $a+b+c+d+e<=10$ ，那么总方案数为： $C_{10+5}^{5}=3003$ 。

2）减去每个超级神器单独超出限制的组合数

对于每个超级神器 $Ti$ 数量为 $Bi$ 件的限制，计算出不满足条件的组合数，并全部减去。还用上例说明：

第1种神器限制数量为1，即不满足条件为 $B1\geq 2$ ，令 $a^{'}=a-2$ ，则问题转化为： $a^{'}+b+c+d+e\leq 10-2$ ，其组合方案为： $C_{8+5}^{5}=1287$ 。同理第2、3、4种神器不满足条件的方案数为： $C_{7+5}^{5}=792,C_{6+5}^{5}=462,C_{5+5}^{5}=252$ ，共2793种。

3）加上同时超过两个限制条件的组合数

对于同时超过两个限制条件的组合共有 $C^{_{4}^{2}}=6$ 种，比如{1,2}神器的组合不满足条件分别为2和3，令 $a^{'}=a-2,b^{'}=b-3$ ，则有 $a^{'}+b^{'}+c+d+e\leq 10-2-3$ ，组合方案数为： $C^{_{10-2-3+5}^{5}}=C^{_{10}^{5}}=252$ ，同理可以计算出另5种{1,3}、{1,4}、{2,3}、{2,4}、{3,4}组合数为： $C^{_{9}^{5}}=126,C^{_{8}^{5}}=56,C^{_{8}^{5}}=56,C^{_{7}^{5}}=21,C^{_{6}^{5}}=6$ ，共517种。

4）减去同时超过三个限制条件的组合数

同时超过三条限制条件的组合，同理有 $C^{_{4}^{4}}=4$ 种，比如{1,2,3}神器的组合不满足条件方程为 $a^{'}+b^{'}+c^{'}+d+e\leq 10-2-3-4$ ,组合方案数为： $C^{_{10-2-3-4+5}^{5}}=C^{_{6}^{5}}=6$ 。{1,2,4}神器的组合数为： $C^{_{10-2-3-5+5}^{5}}=C^{_{5}^{5}}=1$ ，而另两种组合因相减后变负数，则对答案没有贡献，不用解，即方案为0，则共7种。

5）加上同时超过四个限制的组合数

同样，因同时超过四个限制条件的组合，已经超过示例中最多拿10个宝具的总数，即对答案没有贡献，方案为0。

对于T不大于15种的神器均用类似容斥原理进行加减操作，最终结果就是：3003−2793+517−7+0=720种。

具体操作：

预处理阶乘数组及逆元数组

因需要运行 $2^{T}$ (最高是 $2^{15}=32768$ )次计算组合数，建议预处理小于模P下的所有可能的阶乘和逆元数组，以提高计算效率。对于预处理逆元数组，可以使用费马定理+快速幂算法，或用递推公式线性求逆元。建议使用前者，实测效率更高。

费马定理+快速幂算法：

int pow_mod(int a, int b, int P) {// 快速幂求a^b mod pint res = 1;while (b > 0) {if (b & 1) res = (ll)res * a % P;a = (ll)a * a % P;b >>= 1;}return res;
}
void init1(int P) { //预处理阶乘及其逆元1---费马定理+快速幂fact[0] = 1;for (int i = 1; i < P; ++i)fact[i] = (ll)fact[i - 1] * i % P; //阶乘数组invf[P - 1] = pow_mod(fact[P - 1], P - 2, P); //费马定理初始化逆元for (int i = P - 2; i >= 0; --i)invf[i] = (ll)invf[i + 1] * (i + 1) % P; //计算阶乘逆元组
}

递推公式线性求逆元：

void init2(int P) { //预处理阶乘及其逆元2----递推公式线性求逆元fact[0] = invf[0] = invf[1] = 1;for (int i = 2; i < P; i++)invf[i] = (P - (ll)(P / i) * invf[P % i]) % P; //递推公式线性求逆元for (int i = 1; i < P; i++) {fact[i] = (ll)fact[i - 1] * i % P; //阶乘数组invf[i] = (ll)invf[i - 1] * invf[i] % P; //阶乘逆元数组}
}

卢卡斯定理

由于题意给出的模数 $P\leq 10^{5}$ ，小于 $N,M\leq 10^{9}$ ，预处理的阶乘及阶乘逆元数组只能少于 $10^{5}$ ，不能满足求所有组合数，故需要用卢卡斯定理来计算大于模素数P的组合数，即：将M、N展开为P进制下的和形式，再用递归计算出每一位的组合数，最后将所有位的组合数相乘，并对 P 取模求出。

int comb(int n, int k, int P) {//计算组合数 C(n, k) % Pif (k > n) return 0;return (ll)fact[n] * ((ll)invf[k] * invf[n - k] % P) % P;
}
int lucas(int n, int k, int P) { //卢卡斯定理计算n,k大于P的情况if (n < P && k < P) return comb(n, k, P); //小于P直接计算return (ll)lucas(n / P, k / P, P) * comb(n % P, k % P, P) % P; //大于P情况
}

容斥原理求结果

最后通过位掩码循环 $2^{T}$ 次，即枚举所有可能的子集数量。再根据容斥原理求出每一步的组合数的操作结果，即为答案。

附AC代码及截图：

目前在洛谷本题的最优解第一名，总耗时62ms为费马定理+快速幂的预处理方案。

#include <iostream>
#include <vector>
using namespace std;
typedef long long ll;
const int maxn = 100005;
int fact[maxn], invf[maxn];int pow_mod(int a, int b, int P) {// 快速幂求a^b mod pint res = 1;while (b > 0) {if (b & 1) res = (ll)res * a % P;a = (ll)a * a % P;b >>= 1;}return res;
}
void init1(int P) { //预处理阶乘及其逆元1---费马定理+快速幂fact[0] = 1;for (int i = 1; i < P; ++i)fact[i] = (ll)fact[i - 1] * i % P; //阶乘invf[P - 1] = pow_mod(fact[P - 1], P - 2, P); //费马定理初始化逆元for (int i = P - 2; i >= 0; --i)invf[i] = (ll)invf[i + 1] * (i + 1) % P; //阶乘逆元
}
void init2(int P) { //预处理阶乘及其逆元2----递推公式线性求逆元fact[0] = invf[0] = invf[1] = 1;for (int i = 2; i < P; i++)invf[i] = (P - (ll)(P / i) * invf[P % i]) % P; //递推公式线性求逆元for (int i = 1; i < P; i++) {fact[i] = (ll)fact[i - 1] * i % P; //阶乘invf[i] = (ll)invf[i - 1] * invf[i] % P; //阶乘逆元}
}
int comb(int n, int k, int P) {//计算组合数 C(n, k) % Pif (k > n) return 0;return (ll)fact[n] * ((ll)invf[k] * invf[n - k] % P) % P;
}
int lucas(int n, int k, int P) { //卢卡斯定理计算n,k大于P的情况if (n < P && k < P) return comb(n, k, P); //小于P直接计算return (ll)lucas(n / P, k / P, P) * comb(n % P, k % P, P) % P; //大于P情况
}
int Iep(int N, int M, const vector<int> &B, int P) {//容斥原理计算方案数int result = 0;int T = B.size();for (int mask = 0; mask < (1 << T); ++mask) {//位掩码循环,生成所有子集int new_M = M;int sign = 1;for (int i = 0; i < T; ++i) {if (mask & (1 << i)) {//检查mask的二进制表示中第i位是否为1new_M -= B[i] + 1;sign *= -1;}}if (new_M >= 0)result = (result + sign * lucas(N + new_M, new_M, P)) % P;}return (result + P) % P; // 确保结果非负
}int main() {int N, T, M, P;scanf("%d%d%d%d", &N, &T, &M, &P);vector<int> B(T);for (int i = 0; i < T; ++i) scanf("%d", &B[i]);init1(P);//比init2(P)更快！printf("%d\n", Iep(N, M, B, P));return 0;
}

谢谢观赏！

查看全文

http://www.mrgr.cn/news/57277.html