c++ - entera - Rápido n elija k mod p para n grande?

division entera c (2)

Entonces, aquí está cómo puedes resolver tu problema.

Por supuesto que sabes la fórmula:

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

(Ver http://en.wikipedia.org/wiki/Binomial_coefficient#Computing_the_value_of_binomial_coefficients )

Usted sabe cómo calcular el numerador:

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

Ahora, como p es primo, el recíproco de cada entero que es coprime con p está bien definido, es decir, se puede encontrar un ^-1 . Y esto se puede hacer usando el teorema de Fermat a ^p-1 = 1 (mod p) => a * a ^p-2 = 1 (mod p) y entonces a ^-1 = a ^p-2 . Ahora todo lo que necesita hacer es implementar una exponenciación rápida (por ejemplo, utilizando el método binario):

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

Y ahora puedes agregar el denominador a nuestro resultado:

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

Tenga en cuenta que estoy usando long long en todas partes para evitar el desbordamiento de tipos. Por supuesto, no es necesario hacer k exponenciaciones: puede calcular k! (Mod p) y luego dividir solo una vez:

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

EDITAR: según el comentario de @dbaupp si k> = p la k! será igual a 0 módulo p y (k!) ^ - 1 no será definido. Para evitar que primero se compute el grado con el que p está en n * (n-1) ... (n-k + 1) y en k! y compararlos:

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

EDITAR: Hay una optimización más que se puede agregar a la solución anterior: en lugar de calcular el número inverso de cada múltiplo en k !, podemos calcular k! (Mod p) y luego calcular el inverso de ese número. Por lo tanto, tenemos que pagar el logaritmo por la exponenciación solo una vez. Por supuesto, nuevamente tenemos que descartar los divisores p de cada múltiplo. Solo tenemos que cambiar el último ciclo con esto:

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;

Lo que quiero decir con "n grande" es algo en millones. p es primo

He intentado http://apps.topcoder.com/wiki/display/tc/SRM+467 Pero la función parece ser incorrecta (la probé con 144, elijo 6 mod 5 y me da 0 cuando debería darme 2)

Lo he intentado http://online-judge.uva.es/board/viewtopic.php?f=22&t=42690 Pero no lo entiendo completamente

También hice una función recursiva memorizada que usa la lógica (combinaciones (n-1, k-1, p)% p + combinaciones (n-1, k, p)% p) pero me da problemas de desbordamiento de pila porque n es grande

Probé el Teorema de Lucas, pero parece ser lento o inexacto.

Todo lo que estoy tratando de hacer es crear una n rápida / precisa; elija k mod p para n grande. Si alguien pudiera ayudarme a mostrar una buena implementación para esto, estaría muy agradecido. Gracias.

Según lo solicitado, la versión memoial que acierta la pila se desborda para n grande:

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

Para k grande, podemos reducir el trabajo significativamente al explotar dos hechos fundamentales:

Si p es un primo, el exponente de p en la factorización prima de n! está dado por (n - s_p(n)) / (p-1) , donde s_p(n) es la suma de los dígitos de n en la representación base p (entonces para p = 2 , es popcount). Así, el exponente de p en la factorización prima de choose(n,k) es (s_p(k) + s_p(nk) - s_p(n)) / (p-1) , en particular, es cero si y solo si la adición k + (nk) no tiene acarreo cuando se realiza en la base p (el exponente es el número de acarreos).
Teorema de Wilson: p es un primo, si y solo si (p-1)! ≡ (-1) (mod p) (p-1)! ≡ (-1) (mod p) .

El exponente de p en la factorización de n! generalmente se calcula por

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

La verificación de la divisibilidad de choose(n,k) por p no es estrictamente necesaria, pero es razonable tenerla primero, ya que a menudo será el caso, y entonces es menos trabajo:

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

Ahora echemos un vistazo más de cerca a n! . Separamos los números ≤ n en los múltiplos de p y los números coprimos a p . Con

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

Los múltiplos de p contribuyen p^q * q! . Los números coprime p contribuyen al producto de (j*p + k), 1 ≤ k < p para 0 ≤ j < q , y el producto de (q*p + k), 1 ≤ k ≤ r .

Para los números coprime a p solo nos interesará el módulo de contribución p . Cada una de las ejecuciones completas j*p + k, 1 ≤ k < p es congruente con (p-1)! módulo p , por lo que en conjunto producen una contribución de (-1)^q módulo p . La última ejecución (posiblemente) incompleta produce r! módulo p .

Entonces si escribimos

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

obtenemos

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

donde cop(m,r) es el producto de todos los números coprime a p que son ≤ m*p + r .

Hay dos posibilidades, a = b + c y A = B + C , o a = b + c + 1 y A = B + C - p .

En nuestro cálculo, hemos eliminado la segunda posibilidad de antemano, pero eso no es esencial.

En el primer caso, los poderes explícitos de p cancelan, y nos quedamos con

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))

Cualquier poder de p dividir choose(n,k) proviene de choose(a,b) - en nuestro caso, no habrá ninguno, ya que hemos eliminado estos casos antes - y, aunque cop(a,A) / (cop(b,B) * cop(c,C)) no necesita ser un número entero (considere, por ejemplo, choose(19,9) (mod 5) ), cuando se considera la expresión módulo p , cop(m,r) reduce a (-1)^m * r! , entonces, dado que a = b + c , el (-1) cancela y nos quedamos con

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

En el segundo caso, encontramos

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

ya que a = b + c + 1 . El acarreo en el último dígito significa que A < B , entonces módulo p

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

(donde podemos reemplazar la división con una multiplicación por la inversa modular, o verla como una congruencia de números racionales, lo que significa que el numerador es divisible por p ). De todos modos, volvemos a encontrar

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

Ahora podemos recurrir para la parte de choose(a,b) .

Ejemplo:

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)

Ahora la implementación:

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

Para calcular la inversa modular, puede usar el teorema de Fermat (llamado pequeño)

Si p es primo y no es divisible por p , entonces a^(p-1) ≡ 1 (mod p) .

y calcule el inverso como a^(p-2) (mod p) , o utilice un método aplicable a un rango más amplio de argumentos, el algoritmo euclidiano extendido o la expansión continua de fracciones, que le dan la inversa modular para cualquier par de coprime ( positivo) enteros:

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

Como calcular a^(p-2) (mod p) , este es un algoritmo O(log p) , para algunas entradas es significativamente más rápido (en realidad es O(min(log k, log p)) , por lo que para k pequeña y p grande, es considerablemente más rápido), para otros es más lento.

En general, de esta manera necesitamos calcular como máximo O (log_p k) coeficientes binomiales módulo p , donde cada coeficiente binomial necesita como máximo O (p) operaciones, produciendo una complejidad total de operaciones O (p * log_p k). Cuando k es significativamente mayor que p , es mucho mejor que la solución O(k) . Para k <= p , se reduce a la solución O(k) con cierta sobrecarga.