Сколько чисел ниже N взаимно просто с N?

Question

Сколько чисел ниже 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, если вы найдете что-то непонятное, просто спросите в комментариях. Мне интересны алгоритм и математика, идея.

21

algorithm math

Источник

user21384 19 июн '09 в 17:09

4 ответа

Решение

Это функция Эйлера, фи.

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

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

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

30

Источник

user13005 19 июн '09 в 17:11

Вот простая, простая реализация формулы, приведенной на странице википедии, с использованием 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), что представляется вполне приемлемой производительностью.

5

Источник

user95810 19 июн '09 в 18:23

Вот некоторые ссылки на другие обсуждения по этому вопросу, включая некоторые другие языковые реализации:

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

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

1

Источник

user121452 19 июн '09 в 20:13

Другие вопросы по тегам algorithm math

user20713 19 июн '09 в 18:24 2009-06-19 18:24 · Accepted Answer · 2009-06-19 18:24

[Править] Еще одна мысль, которая (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$ ,

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

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

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

Существуют лучшие реализации функции Мёбиуса, и ее можно запоминать по скорости, но за ней должно быть достаточно легко следовать.

Более очевидное вычисление 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'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 с