¿Por qué pow (a, d, n) es mucho más rápido que a ** d% n?

110

Estaba tratando de implementar una prueba de primalidad Miller-Rabin y me sorprendió por qué tardaba tanto (> 20 segundos) en números medianos (~ 7 dígitos). Finalmente encontré que la siguiente línea de código es la fuente del problema:

x = a**d % n

(donde a,, dy nson todos similares, pero desiguales, números medianos, **es el operador de exponenciación y %es el operador de módulo)

Luego intenté reemplazarlo con lo siguiente:

x = pow(a, d, n)

y en comparación es casi instantáneo.

Para el contexto, aquí está la función original:

from random import randint

def primalityTest(n, k):
    if n < 2:
        return False
    if n % 2 == 0:
        return False
    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d >>= 1
    for i in range(k):
        rand = randint(2, n - 2)
        x = rand**d % n         # offending line
        if x == 1 or x == n - 1:
            continue
        for r in range(s):
            toReturn = True
            x = pow(x, 2, n)
            if x == 1:
                return False
            if x == n - 1:
                toReturn = False
                break
        if toReturn:
            return False
    return True

print(primalityTest(2700643,1))

Un ejemplo de cálculo cronometrado:

from timeit import timeit

a = 2505626
d = 1520321
n = 2700643

def testA():
    print(a**d % n)

def testB():
    print(pow(a, d, n))

print("time: %(time)fs" % {"time":timeit("testA()", setup="from __main__ import testA", number=1)})
print("time: %(time)fs" % {"time":timeit("testB()", setup="from __main__ import testB", number=1)})

Salida (ejecutar con PyPy 1.9.0):

2642565
time: 23.785543s
2642565
time: 0.000030s

Salida (ejecutar con Python 3.3.0, 2.7.2 devuelve tiempos muy similares):

2642565
time: 14.426975s
2642565
time: 0.000021s

Y una pregunta relacionada, ¿por qué este cálculo es casi dos veces más rápido cuando se ejecuta con Python 2 o 3 que con PyPy, cuando normalmente PyPy es mucho más rápido ?

Lyallcooper
fuente

Respuestas:

164

Consulte el artículo de Wikipedia sobre exponenciación modular . Básicamente, cuando lo hace a**d % n, en realidad tiene que calcular a**d, lo que podría ser bastante grande. Pero hay formas de computar a**d % nsin tener que computarse a a**dsí mismo, y eso es lo que powhace. El **operador no puede hacer esto porque no puede "ver el futuro" para saber que va a tomar el módulo de inmediato.

BrenBarn
fuente
14
+1 eso es en realidad lo que implica la cadena de documentos>>> print pow.__doc__ pow(x, y[, z]) -> number With two arguments, equivalent to x**y. With three arguments, equivalent to (x**y) % z, but may be more efficient (e.g. for longs).
Hedde van der Heide
6
Dependiendo de su versión de Python, esto solo puede ser cierto bajo ciertas condiciones. IIRC, en 3.xy 2.7, solo puede usar la forma de tres argumentos con tipos integrales (y potencia no negativa), y siempre obtendrá exponenciación modular con el inttipo nativo , pero no necesariamente con otros tipos integrales. Pero en versiones anteriores había reglas sobre cómo encajar en una C long, se permitía la forma de tres argumentos float, etc. (Con suerte, no está usando 2.1 o anterior, y no está usando ningún tipo integral personalizado de módulos C, por lo que ninguno de esto es importante para usted.)
abarnert
13
Según su respuesta, parece que es imposible que un compilador vea la expresión y la optimice, lo cual no es cierto. Es simplemente ocurre que no hay compiladores actuales Python lo hacen.
danielkza
5
@danielkza: Eso es cierto, no quise dar a entender que es teóricamente imposible. Quizás "no mira hacia el futuro" sería más exacto que "no puede ver el futuro". Sin embargo, tenga en cuenta que la optimización podría ser extremadamente difícil o incluso imposible en general. Para constantes operandos que podría ser optimizado, pero en x ** y % n, xpodría ser un objeto que implementa __pow__y, en base a un número aleatorio, devuelve uno de varios objetos diferentes de aplicación __mod__de manera que también dependen de números aleatorios, etc.
BrenBarn
2
@danielkza: Además, las funciones no tienen el mismo dominio: .3 ** .4 % .5es perfectamente legal, pero si el compilador lo transforma en pow(.3, .4, .5)eso generaría un TypeError. El compilador debería poder saber eso a, dy nse garantiza que son valores de un tipo integral (o tal vez solo específicamente de tipo int, porque la transformación no ayuda de otra manera), y dse garantiza que no serán negativos. Eso es algo que un JIT podría hacer, pero un compilador estático para un lenguaje con tipos dinámicos y sin inferencia simplemente no puede.
abarnert
37

BrenBarn respondió a su pregunta principal. Para tu aparte:

¿Por qué es casi el doble de rápido cuando se ejecuta con Python 2 o 3 que con PyPy, cuando normalmente PyPy es mucho más rápido?

Si lee la página de rendimiento de PyPy , este es exactamente el tipo de cosas en las que PyPy no es bueno; de hecho, el primer ejemplo que dan:

Los malos ejemplos incluyen hacer cálculos con largos largos, lo cual se realiza mediante un código de soporte no optimizable.

Teóricamente, convertir una exponenciación enorme seguida de un mod en una exponenciación modular (al menos después de la primera pasada) es una transformación que un JIT podría hacer ... pero no el JIT de PyPy.

Como nota al margen, si necesita hacer cálculos con números enteros grandes, es posible que desee mirar módulos de terceros como gmpy, que a veces puede ser mucho más rápido que la implementación nativa de CPython en algunos casos fuera de los usos principales, y también tiene mucho de funcionalidad adicional que de otro modo tendría que escribir usted mismo, a costa de ser menos conveniente.

abarnert
fuente
2
los largos se arreglaron. pruebe pypy 2.0 beta 1 (no será más rápido que CPython, pero tampoco debería ser más lento). gmpy no tiene una forma de manejar MemoryError :(
fijal
@fijal: Sí, y gmpytambién es más lento en lugar de más rápido en algunos casos, y hace que muchas cosas simples sean menos convenientes. No siempre es la respuesta, pero a veces lo es. Por lo tanto, vale la pena analizarlo si se trata de números enteros grandes y el tipo nativo de Python no parece lo suficientemente rápido.
abarnert
1
y si no le importa si sus números son grandes, su programa se segrega por defecto
fijal
1
Es el factor que hizo que PyPy no usara la biblioteca GMP durante mucho tiempo. Puede que esté bien para usted, pero no para los desarrolladores de Python VM. El malloc puede fallar sin usar mucha RAM, solo coloque un número muy grande allí. El comportamiento de GMP a partir de ese momento no está definido y Python no puede permitirlo.
fijal
1
@fijal: Estoy completamente de acuerdo en que no debería usarse para implementar el tipo integrado de Python. Eso no significa que nunca deba usarse para nada.
abarnert
11

Hay atajos para hacer una exponenciación modular: por ejemplo, puede encontrar a**(2i) mod npara cada ifrom 1to log(d)y multiplicar (mod n) los resultados intermedios que necesita. Una función de exponenciación modular dedicada como 3 argumentos pow()puede aprovechar estos trucos porque sabe que estás haciendo aritmética modular. El analizador de Python no puede reconocer esto dada la expresión simple a**d % n, por lo que realizará el cálculo completo (que llevará mucho más tiempo).

atomicinf
fuente
3

La forma en que x = a**d % nse calcula es subir aa la dpotencia, luego modulo eso con n. En primer lugar, si aes grande, crea un número enorme que luego se trunca. Sin embargo, x = pow(a, d, n)lo más probable es que esté optimizado para que solo nse rastreen los últimos dígitos, que son todo lo que se requiere para calcular el módulo de multiplicación de un número.

Yuushi
fuente
6
"Requiere d multiplicaciones para calcular x ** d" - no es correcto. Puedes hacerlo en multiplicaciones O (log d) (muy amplias). La exponenciación por cuadratura se puede utilizar sin módulo. El gran tamaño de los multiplicandos es lo que lleva la delantera aquí.
John Dvorak
@JanDvorak Es cierto, no estoy seguro de por qué pensé que Python no utilizaría el mismo algoritmo de exponenciación **para pow.
Yuushi
5
No los últimos "n" dígitos ... solo mantiene los cálculos en Z / nZ.
Thomas