快nselectk mod p为大n？

我的意思是“大n”是百万分之一。 p是素数。

我试过http://apps.topcoder.com/wiki/display/tc/SRM+467但function似乎是不正确的（我testing了144select6模5，它给了我0当它应该给我2）

我已经尝试http://online-judge.uva.es/board/viewtopic.php?f=22&t=42690但我完全不明白

我也做了一个使用逻辑（组合（n-1，k-1，p）％p +组合（n-1，k，p）％p）的recursion函数，但是它给了我堆栈溢出的问题， n很大

我已经尝试过卢卡斯定理，但它似乎是缓慢或不准确的。

我所要做的就是创build一个快速/准确的nselectk mod p为大n。 如果有人能够帮助我展示一个很好的实现，我将非常感激。 谢谢。

根据要求，堆栈溢出的memoized版本溢出大n：

std::map<std::pair<long long, long long>, long long> memo; long long combinations(long long n, long long k, long long p){ if (n < k) return 0; if (0 == n) return 0; if (0 == k) return 1; if (n == k) return 1; if (1 == k) return n; map<std::pair<long long, long long>, long long>::iterator it; if((it = memo.find(std::make_pair(n, k))) != memo.end()) { return it->second; } else { long long value = (combinations(n-1, k-1,p)%p + combinations(n-1, k,p)%p)%p; memo.insert(std::make_pair(std::make_pair(n, k), value)); return value; } } 

Answer 1

所以，这里是你如何解决你的问题。

当然你知道这个公式：

 comb(n,k) = n!/(k!*(nk)!) = (n*(n-1)*...(n-k+1))/k! 

（见http://en.wikipedia.org/wiki/Binomial_coefficient#Computing_the_value_of_binomial_coefficients ）

你知道如何计算分子：

 long long res = 1; for (long long i = n; i > n- k; --i) { res = (res * i) % p; } 

现在，因为p是素数，所以与p互质的每个整数的倒数都被很好地定义，即可以find^-1 。 这可以通过使用费马定理a ^p-1 = 1（mod p）=> a * a ^p-2 = 1（mod p）来完成，所以a ^-1 = a ^p-2 。 现在，您只需要执行快速取幂（例如使用二进制方法）：

 long long degree(long long a, long long k, long long p) { long long res = 1; long long cur = a; while (k) { if (k % 2) { res = (res * cur) % p; } k /= 2; cur = (cur * cur) % p; } return res; } 

现在，您可以将分母添加到我们的结果中：

 long long res = 1; for (long long i = 1; i <= k; ++i) { res = (res * degree(i, p- 2)) % p; } 

请注意，我正在使用long long来避免types溢出。 当然，你不需要做k幂运算 – 你可以计算k！（mod p），然后只除以一次：

 long long denom = 1; for (long long i = 1; i <= k; ++i) { denom = (denom * i) % p; } res = (res * degree(denom, p- 2)) % p; 

编辑：根据@ dbaupp的评论，如果K> = P K！ 将等于0模p，（k！）^ – 1将不被定义。 为了避免这种情况，首先计算p在n *（n-1）…（n-k + 1）中的度数， 并比较它们：

 int get_degree(long long n, long long p) { // returns the degree with which p is in n! int degree_num = 0; long long u = p; long long temp = n; while (u <= temp) { degree_num += temp / u; u *= p; } return degree_num; } long long combinations(int n, int k, long long p) { int num_degree = get_degree(n, p) - get_degree(n - k, p); int den_degree = get_degree(k, p); if (num_degree > den_degree) { return 0; } long long res = 1; for (long long i = n; i > n - k; --i) { long long ti = i; while(ti % p == 0) { ti /= p; } res = (res * ti) % p; } for (long long i = 1; i <= k; ++i) { long long ti = i; while(ti % p == 0) { ti /= p; } res = (res * degree(ti, p-2, p)) % p; } return res; } 

编辑：有一个更多的优化，可以添加到上面的解决scheme – 而不是计算在k！每个倍数的倒数，我们可以计算k！（mod p），然后计算该数字的倒数。 因此，我们只需要为求幂支付一次对数。 当然，我们必须抛弃每个倍数的p因数。 我们只需要改变最后一个循环：

 long long denom = 1; for (long long i = 1; i <= k; ++i) { long long ti = i; while(ti % p == 0) { ti /= p; } denom = (denom * ti) % p; } res = (res * degree(denom, p-2, p)) % p; 

Answer 2

对于大k ，我们可以通过利用两个基本事实来显着减less工作：

1. 如果p是素数，则在n!的素数因子分解中p的指数n! (n - s_p(n)) / (p-1) ，其中s_p(n)是基数p表示中的n的数字的总和（因此对于p = 2 ，它是popup数）。 choose(n,k)的素因子化中p的指数为(s_p(k) + s_p(nk) - s_p(n)) / (p-1) ，特别是当且仅当当基数p （指数是进位数）执行时，加法k + (nk)没有进位。

2. 威尔逊定理： p是一个素数，当且仅当(p-1)! ≡ (-1) (mod p) (p-1)! ≡ (-1) (mod p) 。

在n!的因式分解中p的指数 通常是通过计算

 long long factorial_exponent(long long n, long long p) { long long ex = 0; do { n /= p; ex += n; }while(n > 0); return ex; } 

用p来choose(n,k)的可分性检查并不是绝对必要的，但是首先要有这个检查是合理的，因为通常情况是这样的，

 long long choose_mod(long long n, long long k, long long p) { // We deal with the trivial cases first if (k < 0 || n < k) return 0; if (k == 0 || k == n) return 1; // Now check whether choose(n,k) is divisible by p if (factorial_exponent(n) > factorial_exponent(k) + factorial_exponent(nk)) return 0; // If it's not divisible, do the generic work return choose_mod_one(n,k,p); } 

现在让我们仔细看看n! 。 我们把数字≤ n分成p的倍数，数字与p相等。 同

 n = q*p + r, 0 ≤ r < p 

p的倍数有助于p^q * q! 。 与(j*p + k), 1 ≤ k < p的乘积共同构成p的乘积，其中0 ≤ j < q ， 0 ≤ j < q ， (q*p + k), 1 ≤ k ≤ r的乘积(q*p + k), 1 ≤ k ≤ r 。

对于这些数字来说，我们只会对模p的贡献感兴趣。 每个完整运行j*p + k, 1 ≤ k < p与(p-1)! 模p ，所以它们总共产生(-1)^q模p的贡献。 最后（可能）不完整的运行产生r! 模p 。

所以，如果我们写

 n = a*p + A k = b*p + B nk = c*p + C 

我们得到

 choose(n,k) = p^a * a!/ (p^b * b! * p^c * c!) * cop(a,A) / (cop(b,B) * cop(c,C)) 

其中cop(m,r)是所有数的互素与p的乘积，它们≤m ≤ m*p + r 。

有两种可能性： a = b + c和A = B + C ，或a = b + c + 1和A = B + C - p 。

在我们的计算中，我们预先排除了第二种可能性，但这不是必需的。

在第一种情况下， p取消的显式权力，我们留下

 choose(n,k) = a! / (b! * c!) * cop(a,A) / (cop(b,B) * cop(c,C)) = choose(a,b) * cop(a,A) / (cop(b,B) * cop(c,C)) 

choose(n,k)来自choose(a,b) – 在我们的例子中，将不存在，因为我们之前已经消除了这些情况 – 尽pipecop(a,A) / (cop(b,B) * cop(c,C))不必是一个整数（例如考虑choose(19,9) (mod 5) ），当考虑expression式p ， cop(m,r) (-1)^m * r! ，所以，由于a = b + c ， (-1)取消，我们留下

 choose(n,k) ≡ choose(a,b) * choose(A,B) (mod p) 

在第二种情况下，我们发现

 choose(n,k) = choose(a,b) * p * cop(a,A)/ (cop(b,B) * cop(c,C)) 

因为a = b + c + 1 。 最后一位的进位意味着A < B ，所以模p

 p * cop(a,A) / (cop(b,B) * cop(c,C)) ≡ 0 = choose(A,B) 

（我们可以用乘以模的逆来代替除法，或者把它看作有理数的同余，也就是说分子可以被p整除）。 无论如何，我们再次find

 choose(n,k) ≡ choose(a,b) * choose(A,B) (mod p) 

现在我们可以重复choose(a,b)部分。

例：

 choose(144,6) (mod 5) 144 = 28 * 5 + 4 6 = 1 * 5 + 1 choose(144,6) ≡ choose(28,1) * choose(4,1) (mod 5) ≡ choose(3,1) * choose(4,1) (mod 5) ≡ 3 * 4 = 12 ≡ 2 (mod 5) choose(12349,789) ≡ choose(2469,157) * choose(4,4) ≡ choose(493,31) * choose(4,2) * choose(4,4 ≡ choose(98,6) * choose(3,1) * choose(4,2) * choose(4,4) ≡ choose(19,1) * choose(3,1) * choose(3,1) * choose(4,2) * choose(4,4) ≡ 4 * 3 * 3 * 1 * 1 = 36 ≡ 1 (mod 5) 

现在执行：

 // Preconditions: 0 <= k <= n; p > 1 prime long long choose_mod_one(long long n, long long k, long long p) { // For small k, no recursion is necessary if (k < p) return choose_mod_two(n,k,p); long long q_n, r_n, q_k, r_k, choose; q_n = n / p; r_n = n % p; q_k = k / p; r_k = k % p; choose = choose_mod_two(r_n, r_k, p); // If the exponent of p in choose(n,k) isn't determined to be 0 // before the calculation gets serious, short-cut here: /* if (choose == 0) return 0; */ choose *= choose_mod_one(q_n, q_k, p); return choose % p; } // Preconditions: 0 <= k <= min(n,p-1); p > 1 prime long long choose_mod_two(long long n, long long k, long long p) { // reduce n modulo p n %= p; // Trivial checks if (n < k) return 0; if (k == 0 || k == n) return 1; // Now 0 < k < n, save a bit of work if k > n/2 if (k > n/2) k = nk; // calculate numerator and denominator modulo p long long num = n, den = 1; for(n = n-1; k > 1; --n, --k) { num = (num * n) % p; den = (den * k) % p; } // Invert denominator modulo p den = invert_mod(den,p); return (num * den) % p; } 

要计算模逆，可以使用费马（所谓的小）定理

如果p是素数且不能被p整除，则a^(p-1) ≡ 1 (mod p) 。

a^(p-2) (mod p) ，或者使用适用于更宽范围参数的方法，扩展的欧几里德algorithm或连续分数展开，这给出了任何一对副本积极）整数：

 long long invert_mod(long long k, long long m) { if (m == 0) return (k == 1 || k == -1) ? k : 0; if (m < 0) m = -m; k %= m; if (k < 0) k += m; int neg = 1; long long p1 = 1, p2 = 0, k1 = k, m1 = m, q, r, temp; while(k1 > 0) { q = m1 / k1; r = m1 % k1; temp = q*p1 + p2; p2 = p1; p1 = temp; m1 = k1; k1 = r; neg = !neg; } return neg ? m - p2 : p2; } 

就像计算a^(p-2) (mod p) ，这是一个O(log p)algorithm，对于一些input它显着更快（实际上是O(min(log k, log p)) ，所以对于小的k大p ，速度相当快），对于其他人来说，速度更慢。

总的来说，我们需要计算至多O（log_p k）二项式系数mod p ，其中每个二项式系数至多需要O（p）运算，产生O（p * log_p k）运算的总复杂度。 当k明显大于p ，这比O(k)解决scheme好得多。 对于k <= p ，它有一些开销就减less到O(k)解。