algorithm - N以下有多少个数是N的互质数？

Question

简而言之:

假设 a 与 b 互质，如果 GCD(a,b) = 1(其中 GCD 代表 great common divisor )，有多少个 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 的人也应该可以阅读代码，如果您发现任何不清楚的地方，请在评论中提问。我对算法和数学感兴趣，这个想法。

Answer 1

[编辑] 最后一个想法，哪个 (IMO) 足够重要，我会把它放在开头:如果你一次收集一堆 totients，你可以避免很多多余的工作。不要费心从大数字开始寻找它们的较小因子——相反，迭代较小的因子并为较大的数字累积结果。

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 毫秒。

---

Euler totient function 上的维基百科页面有一些不错的数学结果。

计算与 的每个除数互质且小于的数:这有一个简单的*映射来计算 中的整数至 , 所以总和是 .

* 由的第二个定义 trivial

这非常适合 Möbius inversion formula 的应用，一个巧妙的技巧，用于反转这种精确形式的总和。

这很自然地引出了代码

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 function 实现, 它可以被内存以提高速度，但这应该很容易理解。

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 存在更好的实现，但这是为了便于阅读。

---

# 辅助函数

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的第一次尝试

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

# OP 的第二次尝试

# I don't understand the algorithm, and just copying it yields inaccurate results

#莫比乌斯反演

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

#传统公式

def totient3(n):
    return int(reduce(mul, (1 - 1.0 / p for p, exp in factorize(n)), n))

#繁体公式，无除法

def totient4(n):
    return reduce(mul, ((p-1) * p ** (exp-1) for p, exp in factorize(n)), 1)

使用此代码计算桌面上从 1 到 9999 的所有数字的总和，计算 5 次运行的平均值，

• totient1 需要永远
• totient2 耗时10s
• totient3 耗时1.3s
• totient4 耗时1.3s