algorithm - ¿Cuántos números debajo de N son coprimes a N?

math

Question

En breve:

Dado que a es coprime a b si GCD (a, b) = 1 (donde GCD significa gran divisor común ), ¿cuántos enteros positivos debajo de N son coprime a N?

¿Hay alguna manera inteligente?

Cosas no necesarias

Aquí está la manera más tonta:

def count_coprime(N):
counter = 0
for n in xrange(1,N):
if gcd(n,N) == 1:
counter += 1
return counter

Funciona, pero es lento, y tonto. Me gustaría usar un algoritmo inteligente y más rápido. Traté de usar factores primos y divisores de N, pero siempre obtengo algo que no funciona con N. más grande.

Creo que el algoritmo debería poder contarlos sin calcularlos todos, como hace el algoritmo más tonto: P

Editar

Parece que he encontrado un trabajo:

def a_bit_more_clever_counter(N):
result = N - 1
factors = []
for factor, multiplicity in factorGenerator(N):
result -= N/factor - 1
for pf in factors:
if lcm(pf, factor) < N:
result += N/lcm(pf, factor) - 1
factors += [factor]
return result

donde mcm es el mínimo común múltiplo. ¿Alguien tiene uno mejor?

Nota

Estoy usando python, creo que el código debería ser legible incluso para quien no conoce python, si encuentra algo que no esté claro, simplemente pregunte en los comentarios. Me interesa el algoritmo y las matemáticas, la idea.

Answer 1

Aquí hay algunos enlaces a otras discusiones sobre esto, incluyendo algunas implementaciones de otros idiomas:

http://www.velocityreviews.com/forums/t459467-computing-eulers-totient-function.html

http://www.google.com/codesearch?q=Euler%27s+totient&hl=en&btnG=Code

Answer 2

Aquí hay una implementación simple y directa de la fórmula dada en la página de wikipedia, usando gmpy para una fácil factorización (estoy sesgado, pero probablemente quieras gmpy si te importa jugar con cosas enteras divertidas en Python ... ;-):

import gmpy
def prime_factors(x):
prime = gmpy.mpz(2)
x = gmpy.mpz(x)
factors = {}
while x >= prime:
newx, mult = x.remove(prime)
if mult:
factors[prime] = mult
x = newx
prime = prime.next_prime()
return factors
def euler_phi(x):
fac = prime_factors(x)
result = 1
for factor in fac:
result *= (factor-1) * (factor**(fac[factor]-1))
return result

Por ejemplo, en mi modesta estación de trabajo, computar euler_phi (123456789) [para lo cual recibo 82260072] toma 937 microsegundos (con Python 2.5; 897 con 2.4), lo que parece ser un rendimiento bastante razonable.

Answer 3

Esta es la función totient de Euler , phi.

Tiene la emocionante propiedad de ser multiplicativo: si gcd (m, n) = 1, entonces phi (mn) = phi (m) phi (n). Y phi es fácil de calcular para potencias de primos, ya que todo lo que está debajo de ellos es coprime excepto los múltiplos de potencias más pequeñas del mismo primo.

Obviamente, la factorización aún no es un problema trivial, pero incluso las divisiones de prueba sqrt (n) (suficientes para encontrar todos los factores primos) superan las aplicaciones n-1 del algoritmo de Euclid.

Si realiza una memorización, puede reducir el costo promedio de computar una gran cantidad de ellos.

Answer 4

[Editar] Un último pensamiento, que (IMO) es lo suficientemente importante para que lo plantee al principio: si está reuniendo un montón de pacientes a la vez, puede evitar una gran cantidad de trabajo redundante. No se moleste en comenzar con números grandes para encontrar sus factores más pequeños; en su lugar, repita los factores más pequeños y acumule resultados para los números más grandes.

class Totient:
def __init__(self, n):
self.totients = [1 for i in range(n)]
for i in range(2, n):
if self.totients[i] == 1:
for j in range(i, n, i):
self.totients[j] *= i - 1
k = j / i
while k % i == 0:
self.totients[j] *= i
k /= i
def __call__(self, i):
return self.totients[i]
if __name__ == ''__main__'':
from itertools import imap
totient = Totient(10000)
print sum(imap(totient, range(10000)))

Esto toma solo 8ms en mi escritorio.

La página de Wikipedia sobre la función principal de Euler tiene algunos buenos resultados matemáticos.

cuenta los números de coprime y más pequeños que cada divisor de : esto tiene un mapeo trivial * para contar los enteros de a , entonces la suma total es .

* por la segunda definicion de trivial

Esto es perfecto para una aplicación de la fórmula de inversión de Möbius , un truco inteligente para invertir sumas de esta forma exacta.

Esto conduce naturalmente al código.

def totient(n):
if n == 1: return 1
return sum(d * mobius(n / d) for d in range(1, n+1) if n % d == 0)
def mobius(n):
result, i = 1, 2
while n >= i:
if n % i == 0:
n = n / i
if n % i == 0:
return 0
result = -result
i = i + 1
return result

Existen mejores implementaciones de la función de Möbius , y podría ser memorizada para la velocidad, pero esto debería ser bastante fácil de seguir.

El cálculo más obvio de la función totiente es

En otras palabras, factorice completamente el número en primos y exponentes únicos, y haga una simple multiplicación a partir de ahí.

from operator import mul
def totient(n):
return int(reduce(mul, (1 - 1.0 / p for p in prime_factors(n)), n))
def primes_factors(n):
i = 2
while n >= i:
if n % i == 0:
yield i
n = n / i
while n % i == 0:
n = n / i
i = i + 1

Nuevamente, existen mejores implementaciones de prime_factors , pero esto está hecho para una fácil lectura.

# helper functions

from collections import defaultdict
from itertools import count
from operator import mul
def gcd(a, b):
while a != 0: a, b = b % a, a
return b
def lcm(a, b): return a * b / gcd(a, b)
primes_cache, prime_jumps = [], defaultdict(list)
def primes():
prime = 1
for i in count():
if i < len(primes_cache): prime = primes_cache[i]
else:
prime += 1
while prime in prime_jumps:
for skip in prime_jumps[prime]:
prime_jumps[prime + skip] += [skip]
del prime_jumps[prime]
prime += 1
prime_jumps[prime + prime] += [prime]
primes_cache.append(prime)
yield prime
def factorize(n):
for prime in primes():
if prime > n: return
exponent = 0
while n % prime == 0:
exponent, n = exponent + 1, n / prime
if exponent != 0:
yield prime, exponent

# OP''s first attempt

def totient1(n):
counter = 0
for i in xrange(1, n):
if gcd(i, n) == 1:
counter += 1
return counter

# OP''s second attempt

# I don''t understand the algorithm, and just copying it yields inaccurate results

# Möbius inversion

def totient2(n):
if n == 1: return 1
return sum(d * mobius(n / d) for d in xrange(1, n+1) if n % d == 0)
mobius_cache = {}
def mobius(n):
result, stack = 1, [n]
for prime in primes():
if n in mobius_cache:
result = mobius_cache[n]
break
if n % prime == 0:
n /= prime
if n % prime == 0:
result = 0
break
stack.append(n)
if prime > n: break
for n in stack[::-1]:
mobius_cache[n] = result
result = -result
return -result

# traditional formula

def totient3(n):
return int(reduce(mul, (1 - 1.0 / p for p, exp in factorize(n)), n))

# traditional formula, no division

def totient4(n):
return reduce(mul, ((p-1) * p ** (exp-1) for p, exp in factorize(n)), 1)

Usando este código para calcular los totients de todos los números del 1 al 9999 en mi escritorio, con un promedio de más de 5 carreras,

• totient1 toma para siempre
• totient2 toma 10s
• totient3 toma 1.3s
• totient4 toma 1.3s