Rápido n elija k mod p para n grande?

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 stack 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 stack se desborda para n grande:

std::map<std::pair, 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>::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; } } 

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

Por supuesto que sabes la fórmula:

 comb(n,k) = n!/(k!*(nk)!) = (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;