How many numbers below N are coprimes to N?

Solution 1:

[Edit] One last thought, which (IMO) is important enough that I'll put it at the beginning: if you're collecting a bunch of totients at once, you can avoid a lot of redundant work. Don't bother starting from large numbers to find their smaller factors -- instead, iterate over the smaller factors and accumulate results for the larger numbers.

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

This takes just 8ms on my desktop.

The Wikipedia page on the Euler totient function has some nice mathematical results.

counts the numbers coprime to and smaller than each divisor of : this has a trivial* mapping to counting the integers from to , so the sum total is .

* by the second definition of trivial

This is perfect for an application of the Möbius inversion formula, a clever trick for inverting sums of this exact form.

This leads naturally to the code

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

There exist better implementations of the Möbius function, and it could be memoized for speed, but this should be easy enough to follow.

The more obvious computation of the totient function is

In other words, fully factor the number into unique primes and exponents, and do a simple multiplication from there.

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

Again, there exist better implementations of prime_factors, but this is meant for easy reading.

# 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)

Using this code to calculate the totients of all numbers from 1 to 9999 on my desktop, averaging over 5 runs,

• totient1 takes forever
• totient2 takes 10s
• totient3 takes 1.3s
• totient4 takes 1.3s

Solution 2:

This is the Euler totient function, phi.

It has the exciting property of being multiplicative: if gcd(m,n) = 1 then phi(mn) = phi(m)phi(n). And phi is easy to calculate for powers of primes, since everything below them is coprime except for the multiples of smaller powers of the same prime.

Obviously factorization is still not a trivial problem, but even sqrt(n) trial divisions (enough to find all the prime factors) beats the heck out of n-1 applications of Euclid's algorithm.

If you memoize, you can reduce the average cost of computing a whole lot of them.

Solution 3:

Here's a simple, straightforward implementation of the formula given on wikipedia's page, using gmpy for easy factorization (I'm biased, but you probably want gmpy if you care about playing with fun integer stuff in 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

For example, on my modest workstation, computing euler_phi(123456789) [for which I get 82260072] takes 937 microseconds (with Python 2.5; 897 with 2.4), which seems to be quite a reasonable performance.