快速n选择k mod p为大n？

Question

我所说的"大n"是数以百万计的东西.p是素数.

我已经尝试了 http://apps.topcoder.com/wiki/display/tc/SRM+467 但是这个功能似乎是不正确的(我用144选择6 mod 5测试了它,当它应该给我时它给了我0 2)

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

我还做了一个memoized递归函数,它使用逻辑(组合(n-1,k-1,p)%p +组合(n-1,k,p)%p)但是它给了我堆栈溢出问题因为n很大

我尝试过卢卡斯定理,但它似乎要么缓慢还是不准确.

我所要做的就是为大n创建一个快速/准确的n选择k mod p.如果有人能帮我展示一个很好的实现,我将非常感激.谢谢.

根据要求,对于大n的命中堆栈溢出的memoized版本:

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!*(n-k)!) = (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是素数的每一个整数的倒数即互质与对被很好地定义即^-1可以找到.并且这可以使用费马定理来完成,^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;
}

请注意我在任何地方都使用很长时间来避免类型溢出.当然你不需要进行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)和k中的程度!并比较它们:

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;
}

编辑:还有一个优化可以添加到上面的解决方案 - 而不是计算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,我们可以通过利用两个基本事实来显着减少工作:

1. 如果p是素数,则p素数因子分解中的指数n!由(n - s_p(n)) / (p-1)下式给出,其中s_p(n)是n基本p表示中数字的总和(因此p = 2,它是popcount).因此的指数p在质因数分解的choose(n,k)就是(s_p(k) + s_p(n-k) - s_p(n)) / (p-1),特别是,这是零,如果且仅当除了k + (n-k)在碱进行时没有进位p(指数为携带的数量).

2. 威尔逊定理: p只有当时是一个素数(p-1)! ? (-1) (mod p).

p因子分解中的指数n!通常由下式计算

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

对于整除的检查choose(n,k)通过p并非绝对必要,但它是合理的,有第一次,因为它会经常发生的情况,然后它的较少的工作:

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(n-k)) 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!.数字互质p贡献的产品(j*p + k), 1 ? k < p的0 ? j < q,和的产品(q*p + k), 1 ? k ? r.

对于互质数字,p我们只对模数贡献感兴趣p.每个完整的运行j*p + k, 1 ? k < p都是(p-1)!模数一致的p,因此它们总共产生(-1)^q模数的贡献p.最后(可能)不完整的运行产生r!模数p.

所以,如果我们写

n   = a*p + A
k   = b*p + B
n-k = 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*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))

任何p划分的权力choose(n,k)来自choose(a,b)- 在我们的情况下,将没有,因为我们之前已经消除了这些情况 - 并且,虽然cop(a,A) / (cop(b,B) * cop(c,C))不需要是整数(例如考虑choose(19,9) (mod 5)),但在考虑表达式模数时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).无论如何,我们再次找到

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 = n-k;
    // 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是素数a而不能被整除p,那么a^(p-1) ? 1 (mod p).

并计算逆a^(p-2) (mod p)α,或使用适用于更广泛的参数的方法,扩展的欧几里德算法或连续的分数扩展,它为任何一对互质(正)整数提供模块化逆:

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)算法,对于某些输入它显着更快(实际上O(min(log k, log p)),对于小型k和大型p,它要快得多),而对于其他输入则更慢.

总的来说,这种方式我们需要计算最多O(log_p k)二项式系数模数p,其中每个二项式系数最多需要O(p)运算,产生O(p*log_p k)运算的总复杂度.当k显着大于时p,这比O(k)解决方案好得多.因为k <= p,它减少O(k)了一些开销的解决方案.