[Редактировать] Еще одна мысль, которая (IMO) достаточно важна, и я поставлю ее в начале: если вы собираете кучу тентов сразу, вы можете избежать многих избыточной работы. Не пытайтесь начинать с больших чисел, чтобы найти их меньшие факторы - вместо этого перебирайте меньшие факторы и накапливайте результаты для больших чисел.
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)))
Это займет всего 8 мс на моем рабочем столе.
На странице Википедии в totient function есть несколько хороших математических результатов.
считает числа, взаимно простые и меньшие, чем каждый делитель 
: это имеет тривиальное * отображение для подсчета целых чисел от 
до , поэтому общая сумма равна .
* по второму определению тривиально
Это идеально подходит для применения формулы обращения Мёбиуса , умного трюка для инвертирования сумм этой точной формы.
Это естественно приводит к коду
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
Существуют более совершенные реализации функции Мёбиуса , и ее можно запоминать по скорости, но это должно быть достаточно просто для отслеживания.
Более очевидное вычисление totient-функции:
Другими словами, полностью разделить число на уникальные простые числа и показатели степени и выполнить простое умножение оттуда.
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
Опять же, существуют лучшие реализации prime_factors, но это предназначено для легкого чтения.
# 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)
Использование этого кода для вычисления числа всех чисел от 1 до 9999 на моем рабочем столе, в среднем за 5 прогонов,
totient1 длится вечно totient2 занимает 10 с totient3 занимает 1,3 с totient4 занимает 1,3 с