algorithm - residuales - ecuaciones de congruencia ejercicios resueltos

Manera rápida de calcular n! mod m donde m es primo

Question

Tenía curiosidad de si había una buena manera de hacer esto. Mi código actual es algo así como:

def factorialMod(n, modulus):
ans=1
for i in range(1,n+1):
ans = ans * i % modulus
return ans % modulus

¡Pero parece bastante lento!

¡Tampoco puedo calcular n! y luego aplicar el módulo primo porque a veces n es tan grande que n! simplemente no es factible calcular explícitamente.

También me encontré con http://en.wikipedia.org/wiki/Stirling%27s_approximation y me pregunto si esto se puede usar aquí de alguna manera.

O bien, ¿cómo podría crear una función recursiva y memorable en C ++?

Answer 1

n puede ser arbitrariamente grande

Bueno, n no puede ser arbitrariamente grande, si n >= m , entonces n! ≡ 0 (mod m) n! ≡ 0 (mod m) (porque m es uno de los factores, por la definición de factorial) .

Asumiendo n << m y necesitas un valor exacto , tu algoritmo no puede ser más rápido, que yo sepa. Sin embargo, si n > m/2 , puede usar la siguiente identidad ( el teorema de Wilson - ¡Gracias @ Daniel Fischer!)

para limitar el número de multiplicaciones a aproximadamente mn

(m-1)! ≡ -1 (mod m) 1 * 2 * 3 * ... * (n-1) * n * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m) n! * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m) n! ≡ -[(n+1) * ... * (m-2) * (m-1)]^-1 (mod m)

Esto nos da una forma simple de calcular n! (mod m) n! (mod m) en mn-1 multiplicaciones, más una inversa modular :

def factorialMod(n, modulus): ans=1 if n <= modulus//2: #calculate the factorial normally (right argument of range() is exclusive) for i in range(1,n+1): ans = (ans * i) % modulus else: #Fancypants method for large n for i in range(n+1,modulus): ans = (ans * i) % modulus ans = modinv(ans, modulus) ans = -1*ans + modulus return ans % modulus

Podemos reformular la ecuación anterior de otra manera, que puede o no puede realizarse un poco más rápido. Usando la siguiente identidad:

podemos reformular la ecuación como

n! ≡ -[(n+1) * ... * (m-2) * (m-1)]^-1 (mod m) n! ≡ -[(n+1-m) * ... * (m-2-m) * (m-1-m)]^-1 (mod m) (reverse order of terms) n! ≡ -[(-1) * (-2) * ... * -(m-n-2) * -(m-n-1)]^-1 (mod m) n! ≡ -[(1) * (2) * ... * (m-n-2) * (m-n-1) * (-1)^(m-n-1)]^-1 (mod m) n! ≡ [(m-n-1)!]^-1 * (-1)^(m-n) (mod m)

Esto se puede escribir en Python de la siguiente manera:

def factorialMod(n, modulus): ans=1 if n <= modulus//2: #calculate the factorial normally (right argument of range() is exclusive) for i in range(1,n+1): ans = (ans * i) % modulus else: #Fancypants method for large n for i in range(1,modulus-n): ans = (ans * i) % modulus ans = modinv(ans, modulus) #Since m is an odd-prime, (-1)^(m-n) = -1 if n is even, +1 if n is odd if n % 2 == 0: ans = -1*ans + modulus return ans % modulus

Si no necesita un valor exacto , la vida se vuelve un poco más fácil: puede usar http://en.wikipedia.org/wiki/Stirling%27s_approximation para calcular un valor aproximado en el tiempo O(log n) (usando la exponenciación cuadrando ) .

Finalmente, debo mencionar que si esto es crítico en el tiempo y está usando Python, intente cambiar a C ++. Por experiencia personal, debe esperar un aumento de velocidad de orden de magnitud o más, simplemente porque este es exactamente el tipo de lazo apretado de CPU que el código compilado de forma nativa sobresale en (también, por alguna razón, GMP parece mucho más afinado que el Bignum de Python) .

Answer 2

¡norte! mod m se puede calcular en O (n ^{1/2 + ε} ) operaciones en lugar de la ingenua O (n). Esto requiere el uso de la multiplicación de polinomios de FFT, y solo vale la pena para n muy grande, por ejemplo, n> 10 ⁴ .

Aquí se puede ver un resumen del algoritmo y algunos tiempos: http://fredrikj.net/blog/2012/03/factorials-mod-n-and-wilsons-theorem/

Answer 3

Ampliando mi comentario a una respuesta:

Sí, hay formas más eficientes de hacer esto. Pero son extremadamente desordenados.

Entonces, a menos que realmente necesite ese rendimiento extra, no sugiero que intente implementarlos.

La clave es observar que el módulo (que es esencialmente una división) va a ser la operación de cuello de botella. Afortunadamente, hay algunos algoritmos muy rápidos que le permiten realizar módulos sobre el mismo número muchas veces.

• División por enteros invariantes utilizando multiplicación
• Reducción de Montgomery

Estos métodos son rápidos porque esencialmente eliminan el módulo.

Esos métodos solo deberían darle una aceleración moderada. Para ser realmente eficiente, es posible que deba desenrollar el ciclo para permitir un mejor IPC:

Algo como esto:

ans0 = 1
ans1 = 1
for i in range(1,(n+1) / 2):
ans0 = ans0 * (2*i + 0) % modulus
ans1 = ans1 * (2*i + 1) % modulus
return ans0 * ans1 % modulus

pero teniendo en cuenta un número impar de iteraciones y combinándolo con uno de los métodos a los que he vinculado anteriormente.

Algunos pueden argumentar que el desenrollado de bucles debe dejarse al compilador. Voy a contraargumentar que los compiladores actualmente no son lo suficientemente inteligentes como para desenrollar este ciclo en particular. Eche un vistazo más de cerca y verá por qué.

Tenga en cuenta que, aunque mi respuesta es independiente del idioma, está destinada principalmente para C o C ++.

Answer 4

Ampliando mi comentario, esto toma alrededor del 50% del tiempo para todo n en [100, 100007] donde m = (117 | 1117):

Function facmod(n As Integer, m As Integer) As Integer
Dim f As Integer = 1
For i As Integer = 2 To n
f = f * i
If f > m Then
f = f Mod m
End If
Next
Return f
End Function

Answer 5

Encontré esta siguiente función en quora:
Con f (n, m) = n! mod m;

function f(n,m:int64):int64;
begin
if n = 1 then f:= 1
else f:= ((n mod m)*(f(n-1,m) mod m)) mod m;
end;

Probablemente supere el uso de un ciclo que consume mucho tiempo y la multiplicación de un gran número almacenado en una cadena. Además, es aplicable a cualquier número entero m.
El enlace donde encontré esta función: https://www.quora.com/How-do-you-calculate-n-mod-m-where-n-is-in-the-1000s-and-m-is-a-very-large-prime-number-eg-n-1000-m-10-9+7

Answer 6

Este algoritmo tiene O(n log log(n)) complejidad de tiempo de preprocesamiento (debido a tamiz) y O (r (n)) complejidad de espacio donde r(n) es la función de conteo principal . Después del preprocesamiento, preguntando x! donde x <= N es O(r(x) * log(x)) .

Desafortunadamente, esto se vuelve inviable para 10 ¹⁰ ya que r (1e9) es 455,052,511 que requieren casi 2 GB de memoria para almacenar los enteros principales. Para 10 ⁹ , sin embargo, solo necesitábamos alrededor de 300 MB de memoria.

Supongamos que mod es 1e9 + 7 para evitar el desbordamiento en C ++, pero puede ser cualquier cosa.

Una forma de calcular N! mod m rápido

Para calcular N! rápido, podemos descomponer N en sus factores primos. Se puede observar que si p es un factor primordial de N !, p debe ser <= N desde N! es un producto de los números enteros de 1 a N. Dados estos hechos, podemos calcular N! usando solo números primos entre 2 a N (inclusive). Según Wikipedia , hay 50,847,534 primos (~ 5 * 10 ^ 7) primos menores o iguales a 10 ⁹ . Esta sería la fórmula para usar:

v _{p (N!)} es el número más grande tal que p ^{v _p (N!)} divide N !. La fórmula de Legendre calcula v _p (N!) En el tiempo O (log (N)). Aquí está la fórmula de Wikipedia.

int factorial(int n) {
int ans = 1;
for (int i = 0; i < primes.size() && primes[i] <= n; i++) {
int e = 0;
for (int j = n; j > 0; ) {
j /= primes[i];
e += j;
}
ans = ((u64) ans * fast_power(primes[i], e)) % mod;
}
return ans;
}

fast_power(x, y) devuelve x ^y % mod en el tiempo O (log (y)) usando exponenciación por cuadratura .

Generando Primes

Estoy generando números primos utilizando Sieve of Eratosthenes . En mi implementación de tamiz, estoy usando bitset para ahorrar memoria. Mi implementación tomó 16-18 segundos para generar números primos hasta 10 ^ 9 cuando compilé con la bandera -O3 con un procesador i3 de 2da generación . Estoy seguro de que hay mejores algoritmos / implementaciones para generar primos.

const int LIMIT = 1e9;
bitset<LIMIT+1> sieve;
vector<int> primes;
void go_sieve() {
primes.push_back(2);
int i = 3;
for (int i = 3; i * i <= LIMIT; i+=2) {
if (!sieve.test(i)) {
primes.push_back(i);
for (int j = i * i; j <= LIMIT; j += i) {
sieve.set(j);
}
}
}
while (i <= LIMIT) {
if (!sieve.test(i)) primes.push_back(i);
i+=2;
}
}

Consumo de memoria

Según Wikipedia , hay 50,847,534 primos que son menores o iguales a 10 ^ 9. Como un entero de 32 bits tiene 4 bytes, el vector principal necesita 203,39 MB (50,847,534 * 4 bytes). El conjunto de bits del tamiz necesita (125 MB) 10 ^ 9 bits.

Actuación

¡Pude calcular 1,000,000,000! en 1.17 segundos después del procesamiento.

Nota: Hice algoritmo con -O3 bandera habilitada.

Answer 7

Si n = (m - 1) para m primero, entonces por http://en.wikipedia.org/wiki/Wilson ''s_theorem n! mod m = (m - 1)

También como ya se ha señalado n! mod m = 0 si n> m

Answer 8

Si queremos calcular M = a*(a+1) * ... * (b-1) * b (mod p) , podemos usar el siguiente enfoque, si suponemos que podemos sumar, restar y multiplicar rápido ( mod p), y obtener una complejidad de tiempo de ejecución de O( sqrt(ba) * polylog(ba) ) .

Para simplificar, supongamos (b-a+1) = k^2 , es un cuadrado. Ahora, podemos dividir nuestro producto en k partes, es decir, M = [a*..*(a+k-1)] *...* [(b-k+1)*..*b] . Cada uno de los factores en este producto es de la forma p(x)=x*..*(x+k-1) , para x apropiado.

Mediante el uso de un algoritmo de multiplicación rápida de polinomios, como el algoritmo Schönhage-Strassen , de manera dividida y vencida, se pueden encontrar los coeficientes del polinomio p(x) in O( k * polylog(k) ) . Ahora, aparentemente hay un algoritmo para sustituir k puntos en el mismo polinomio de grado-k en O( k * polylog(k) ) , lo que significa que podemos calcular p(a), p(a+k), ..., p(b-k+1) rápido.

Este algoritmo de sustituir muchos puntos en un polinomio se describe en el libro "Números primos" por C. Pomerance y R. Crandall. Eventualmente, cuando tenga estos valores k , puede multiplicarlos en O(k) y obtener el valor deseado.

Tenga en cuenta que todas nuestras operaciones fueron tomadas (mod p) . El tiempo exacto de ejecución es O(sqrt(ba) * log(ba)^2 * log(log(ba))) .

Answer 9

Suponiendo que el operador "mod" de tu plataforma elegida sea lo suficientemente rápido, ¡estás limitado principalmente por la velocidad a la que puedes calcular n! y el espacio que tienes disponible para calcularlo en

Entonces, es esencialmente una operación de 2 pasos:

1. Calcula n! (hay muchos algoritmos rápidos, así que no repetiré ninguno aquí)
2. Toma el mod del resultado

No hay necesidad de complejizar las cosas, especialmente si la velocidad es el componente crítico. En general, haga tan pocas operaciones dentro del ciclo como pueda.

Si necesita calcular n! mod m n! mod m repetidamente, entonces tal vez quiera memorizar los valores que salen de la función haciendo los cálculos. Como siempre, es el intercambio clásico de espacio / tiempo, pero las tablas de búsqueda son muy rápidas.

Por último, puede combinar la memorización con la recursión (y trampolines también, si es necesario) para hacer las cosas realmente rápido.