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

division entera c

Question

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

Answer 1

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;

Answer 2

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

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

2. 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.