1

我正在尝试解决黑客等级的编码挑战，该挑战需要计算二项式系数模素数，即

nchoosek(n, k, p)

我正在使用此答案中的代码，该代码适用于前三组输入，但在 4 日开始失败。我在调试器中单步执行并确定问题出现在以下情况：

n % p == 0 || k % p == 0

我只需要知道如何修改我当前的解决方案来处理 n % p == 0 或 k % p == 0 的特定情况。我在堆栈交换中找到的答案似乎都没有解决这个特定情况。这是我的代码：

#include <iostream>
#include <fstream>

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

unsigned long long ModularMultiply(unsigned long long a, unsigned long long b, unsigned long p) {
  unsigned long long a1 = (a >> 21), a2 = a & ((1ull << 21) - 1);
  unsigned long long temp = (a1 * b) % p;   // doesn't overflow under the assumptions
  temp = (temp << 21) % p;                  // this neither
  temp += (a2 * b) % p;                       // nor this
  return temp % p;
}

unsigned long long ModularInverse(unsigned long long k, unsigned 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;
    unsigned 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;
}

// Preconditions: 0 <= k <= min(n,p-1); p > 1 prime
unsigned long long ChooseModTwo(unsigned long long n, unsigned long long k, unsigned 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 = n-k;
    }
    // calculate numerator and denominator modulo p
    unsigned long long num = n, den = 1;
    for(n = n-1; k > 1; --n, --k)
    {
        num = ModularMultiply(num, n, p);
        den = ModularMultiply(den, k, p);
    }

    den = ModularInverse(den,p);
    return ModularMultiply(num, den, p);
}

// Preconditions: 0 <= k <= n; p > 1 prime
long long ChooseModOne(long long n, long long k, const unsigned long p)
{
    // For small k, no recursion is necessary
    if (k < p) return ChooseModTwo(n,k,p);
    unsigned 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 = ChooseModTwo(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;
    return ModularMultiply(choose, ChooseModOne(q_n, q_k, p), p);

}

unsigned long long ModularBinomialCoefficient(unsigned long long n, unsigned long long k, const unsigned 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 (FactorialExponent(n, p) > FactorialExponent(k, p) + FactorialExponent(n - k, p)) return 0;
    // If it's not divisible, do the generic work
    return ChooseModOne(n, k, p);
}

int main() {

  //std::ifstream fin ("input03.txt");
  std::ifstream fin ("test.in");

  int kMod = 1000003;
  int T;
  fin >> T;
  int N = T;
  //std::cin >> T;

  unsigned long long n, k;
  unsigned long long a, b;
  int result[N];
  int index = 0;

  while (T--) {

    fin >> n >> k;

    a = ModularBinomialCoefficient(n - 3, k, kMod);
    b = ModularBinomialCoefficient(n + k, n - 1, kMod);

    // (1 / (n + k) * nCk(n - 3, k) * nCk(n + k, n - 1)) % 1000003
    unsigned long long x = ModularMultiply(a, b, kMod);
    unsigned long long y = ModularMultiply(x, ModularInverse((n + k), kMod), kMod);

    result[index] = y;
    index++;
  }

  for(int i = 0; i < N; i++) {
    std::cout << result[i] << "\n";
  }

  return 0;
}

Input:
6
90 13
65434244 16341234
23424244 12341234
424175 341198
7452123  23472
56000168 16000048

Output:
815483
715724
92308
903465
241972
0 <-- Incorrect, should be: 803478

约束：
4 <= N <= 10^9
1 <= K <= N

4

2 回答 2

2

您可以使用卢卡斯定理将问题简化为 ceil(log_P(N)) 子问题k, n < p：写入n = n_m * p^m + ... + n_0和k = k_m * p^m + ... + k_0在基数p（n_i, k_i < p是数字），然后我们有

C(n,k) = PROD(i = 0 to m, C(n_i, k_i))  (mod p)

子问题很容易解决，因为每个因子k!都有一个反模 p。p << n你会得到一个运行时复杂度为 O(p log(n)) 的算法，如果我理解正确的话，它比 Ivaylo 的代码更好。

int powmod(int x, int e, int p) {
    if (e == 0) return 1;
    if (e & 1) return (long long)x * powmod(x, e - 1, p) % p;
    long long rt = powmod(x, e / 2, p);
    return rt * rt % p;
}

int binom_coeff_mod_prime(int n, int k, int p) {
    long long res = 1;
    while (n || k) {
        int N = n % p, K = k % p;
        for (int i = N - K + 1; i <= N; ++i)
            res = res * i % p;
        for (int i = 1; i <= K; ++i)
            res = res * powmod(i, p - 2, p) % p;
        n /= p;
        k /= p;
    }
    return res;
}
于 2014-05-29T18:59:35.663 回答
0

I suggest you use factorization to compute the number of combinations without division. I've got code for doing so here, originally inspired by Fast computation of multi-category number of combinations (I still would like to post a proper answer to that, if some kind souls would reopen it).

My code stores the result as a table of factors, doing the modular multiplication to expand the result should be quite straightforward.

Probably not practical for n in the range of 10**9, though, since the sieve will be quite massive and take a while to construct.

于 2014-05-29T18:42:46.560 回答