Быстрое n выбирает k mod p для больших n?

Ответ 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 является простым, обратное целое целое, совпадающее с p, хорошо определено, т.е. можно найти a ^-1. И это можно сделать, используя теорему Ферма a ^p-1= 1 (mod p) = > a * a ^p-2= 1 (mod p), и поэтому a ¹= а ^п-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 везде, чтобы избежать переполнения типов. Конечно, вам не нужно делать k exponentiations - вы можете вычислить 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;

EDIT: согласно комментарию @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;
}

EDIT: есть еще одна оптимизация, которая может быть добавлена ​​к решению выше - вместо вычисления обратного числа каждого кратного в 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;

Ответ 2

При больших k мы можем значительно сократить работу, используя два фундаментальных факта:

• Если 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 (показатель степени равен количество переносов).

• Теорема Вильсона: 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 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, поэтому modulo 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), или использовать метод, применимый к более широкому диапазону аргументов, расширенный алгоритм Евклида или продолженное дробное расширение, которое дает вам модульный обратный для любой пары взаимно простых (положительных) целых чисел:/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) с некоторыми накладными расходами.