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

En breve:

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

¿Hay alguna manera inteligente?


No es necesario

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

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

Editar

Parece que encontré uno que funciona:

 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 lcm es el mínimo común múltiplo. ¿Alguien tiene uno mejor?

Nota

Estoy usando Python, creo que el código debe ser legible incluso para quien no conoce Python, si encuentra algo que no está claro solo pregunte en los comentarios. Estoy interesado en el algoritmo y las matemáticas, la idea.

[Editar] Un último pensamiento, que (IMO) es lo suficientemente importante como para ponerlo al principio: si está reuniendo un montón de totients 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, itere sobre los factores más pequeños y acumule los 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 totient de Euler tiene algunos buenos resultados matemáticos.

\ sum_ {d \ mid n} \ varphi (d) cuenta los números coprime y más pequeños que cada divisor de norte : esto tiene un mapeo trivial * para contar los enteros de 1 a norte , entonces la sum total es norte .

* por la segunda definición de trivial

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

\ varphi (n) = \ sum_ {d \ mid n} d \ cdot \ mu \ left (\ frac nd \ right)

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 memorada para la velocidad, pero esto debería ser lo suficientemente fácil de seguir.

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

\ varphi \ left (p_1 ^ {k_1} \ dots p_r ^ {k_r} \ right) = (p_1-1) p_1 ^ {k_1-1} \ dots (p_r-1) p_r ^ {k_r-1} p_1 ^ { k_1} \ dots p_r ^ {k_r} \ prod_ {i = 1} ^ r \ left (1- \ frac1 {p_r} \ right)

En otras palabras, factoriza completamente el número en primos y exponentes únicos, y haz una simple multiplicación desde allí.

 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 es para facilitar la 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, promediando más de 5 carreras,

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

Esta es la función de totler 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 los poderes de primos, ya que todo lo que está debajo de ellos es coprime, excepto por los múltiplos de potencias más pequeñas del mismo primo.

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

Si memoriza, puede reducir el costo promedio de computar muchos de ellos.

Aquí hay una implementación simple y directa de la fórmula dada en la página de wikipedia, usando gmpy para facilitar la factorización (soy parcial, pero probablemente quieras gmpy si te interesa jugar con enteros divertidos 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, calcular euler_phi (123456789) [para el que obtengo 82260072] lleva 937 microsegundos (con Python 2.5; 897 con 2.4), lo que parece ser un rendimiento bastante razonable.

Aquí hay algunos enlaces a otras discusiones sobre esto, incluyendo algunas implementaciones en 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