代码之家 › 专栏 › 技术社区 › Andrea Ambu

n以下有多少个数是n的同模数?

math algorithm

Andrea Ambu · 技术社区 · 16 年前

简而言之:

鉴于一互素对乙如果 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的人来说,代码也应该是可读的,如果您发现任何不清楚的地方,只需在注释中询问即可。我对算法、数学和思想很感兴趣。

4 回复 | 直到 16 年前

Community CDub 8 年前

[编辑] 最后一个想法(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毫秒。

维基百科网页 Euler totient function 有一些很好的数学结果。

$\sum_{d\mid n}\varphi(d)$ 对与每个除数互质且小于的数进行计数。 $\!n$ :这有一个简单的*映射,用于计算 $\!1$ 到 $\!n$ ,所以总数是 $\!n$ .

*根据第二个定义 trivial

这对于 MÃ¶bius inversion formula 这是一个巧妙的技巧,可以将这种精确形式的和求反。

$\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 function 它可以被记忆为速度,但这应该很容易理解。

对tottent函数的更明显的计算是

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

Steve Jessop 16 年前

这就是 Euler totient function ,φ。

它具有乘法的激励性质:如果gcd(m,n)=1,那么phi(m n)=phi(m)phi(n)。phi很容易计算素数的幂,因为除同一素数较小幂的倍数外,下面的一切都是互质的。

显然,因子分解仍然不是一个微不足道的问题,但即使是sqrt(n)试验分区(足以找到所有的素因子)也比欧几里得算法的n-1应用更为困难。

如果你记住了,你就可以降低计算它们的平均成本。

Alex Martelli 16 年前

下面是一个简单、简单的wikipedia页面上给出的公式的实现,使用gmpy进行简单的因子分解(我有偏见,但如果你想在python中玩有趣的整数,你可能需要gmpy…;-):

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),这似乎是相当合理的性能。

iwanttoprogram 16 年前

以下是有关此问题的其他讨论的链接-包括一些其他语言实现:

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

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