Подтвердить что ты не робот

Сколько чисел ниже N являются взаимными ошибками до N?

Короче:

Учитывая, что a является общим для b, если GCD (a, b) = 1 (где GCD означает отличный общий делитель), сколько положительных целых чисел ниже N взаимно просты в N?

Есть ли умный способ?

Не требуется материал

Вот самый тупой способ:

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

Он работает, но он медленный и немой. Я бы хотел использовать умный и быстрый алгоритм. Я попытался использовать простые множители и делители N, но всегда получаю то, что не работает с большим N.

Я думаю, что алгоритм должен иметь возможность подсчитывать их, не вычисляя их, как самый тупой алгоритм: P

Изменить

Кажется, я нашел рабочий:

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

где lcm является наименее общим кратным. Кто-нибудь есть лучший?

Примечание

Я использую python, я думаю, что код должен быть доступен для чтения даже для тех, кто не знает python, если вы найдете что-то, что неясно, просто спросите в комментариях. Меня интересует алгоритм и математика, идея.

4b9b3361

ОТВЕТЫ

Ответ 1

[Edit]. Последняя мысль, которая (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 мс на моем рабочем столе.

На странице Википедии функция тотализатора Эйлера имеет хорошие математические результаты.

\sum_{d\mid n}\varphi(d) подсчитывает числа, равные и меньше каждого делителя n: это имеет тривиальное * отображение для подсчета целых чисел от 1 до n, поэтому сумма равна n.

* вторым определением тривиальным

Это идеально подходит для применения формулы обращения Möbius, умного трюка для инверсии сумм этой точной формы.

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

Это естественно приводит к коду

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

Существуют более реалистичные реализации функции Möbius, и его можно было бы запомнить для скорости, но этого должно быть достаточно просто, чтобы следовать.

Более очевидным вычислением функции totient является

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

Другими словами, полностью умножить число на единственные простые числа и показатели и сделать там простое умножение. 
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 first attempt

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

# OP 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.3s

Ответ 2

Это функция Euler totient, phi.

Он обладает захватывающим свойством быть мультипликативным: если gcd (m, n) = 1, то phi (mn) = phi (m) phi (n). И phi легко вычислить для степеней простых чисел, так как все под ними взаимно просты, за исключением кратных меньших степеней одного и того же простого числа.

Очевидно, что факторизация все еще не является тривиальной проблемой, но даже sqrt (n) пробных делений (достаточно, чтобы найти все основные факторы) превосходит n-1 приложений алгоритма Евклида.

Если вы помните, вы можете уменьшить среднюю стоимость вычисления их большого количества.

Ответ 3

Здесь простая, простая реализация формулы, приведенной на странице wikipedia, с использованием gmpy для легкой факторизации (я предвзятый, но вы, вероятно, хотите gmpy, если вам интересно играть с интересным целым материалом в 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

Например, на моей скромной рабочей станции вычисление euler_phi (123456789) [для которого я получаю 82260072] занимает 937 микросекунд (с Python 2.5; 897 с 2.4), что, по-видимому, вполне разумно.