N 以下有多少个数是 N 的互质数? [英] How many numbers below N are coprimes to N?

问题描述

简而言之:

假设 a 与 b 互质，如果 GCD(a,b) = 1(其中 GCD 代表

这非常适合应用存在更好的实现，它可能是为了速度而记忆，但这应该很容易理解.

totient 函数更明显的计算是

换句话说，将数字完全分解为唯一的素数和指数，然后从那里做一个简单的乘法.

from 操作符 import muldef totient(n):return int(reduce(mul, (1 - 1.0/p for p in prime_factors(n)), n))def primes_factors(n):我 = 2当 n >= i 时:如果 n % i == 0:产量我n = n/i而 n % i == 0:n = n/i我 = 我 + 1

同样，存在更好的 prime_factors 实现，但这是为了便于阅读.

<小时>

# 辅助函数

from collections import defaultdict从 itertools 导入计数从操作符导入 muldef gcd(a, b):而 a != 0: a, b = b % a, a返回 bdef lcm(a, b): 返回 a * b/gcd(a, b)primes_cache, prime_jumps = [], defaultdict(list)定义素数():素数 = 1对于计数()中的我:如果我＜len(primes_cache): 素数 = primes_cache[i]别的:素数 += 1而prime_jumps中的素数:在 prime_jumps[prime] 中跳过:prime_jumps[素数 + 跳过] += [跳过]del prime_jumps[素数]素数 += 1prime_jumps[素数 + 素数] += [素数]primes_cache.append(素数)产量素数定义分解(n):对于素数()中的素数:如果素数 >n:返回指数 = 0而 n % 素数 == 0:指数，n = 指数 + 1，n/素数如果指数 != 0:产量素数，指数

# OP 的第一次尝试

def totient1(n):计数器 = 0对于 xrange(1, n) 中的 i:如果 gcd(i, n) == 1:计数器 += 1退货柜台

# OP 的第二次尝试

#我不懂算法，只是复制它会产生不准确的结果

# 莫比乌斯反演

def totient2(n):如果 n == 1:返回 1返回 sum(d * mobius(n/d) for d in xrange(1, n+1) if n % d == 0)mobius_cache = {}def mobius(n):结果，堆栈 = 1，[n]对于素数()中的素数:如果在 mobius_cache 中为 n:结果 = mobius_cache[n]休息如果 n % 素数 == 0:n/= 素数如果 n % 素数 == 0:结果 = 0休息stack.append(n)如果素数 >n:中断对于堆栈中的 n[::-1]:mobius_cache[n] = 结果结果 = - 结果返回结果

# 传统公式

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 需要 10 秒
• totient3 需要 1.3 秒
• totient4 需要 1.3 秒

In short:

Given that a is coprime to b if GCD(a,b) = 1 (where GCD stands for great common divisor), how many positive integers below N are coprime to N?

Is there a clever way?

---

Not necessary stuff

Here is the dumbest way:

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

It works, but it is slow, and dumb. I'd like to use a clever and faster algorithm. I tried to use prime factors and divisors of N but I always get something that doesn't work with larger N.

I think the algorithm should be able to count them without calculating all of them like the dumbest algorithm does :P

Edit

It seems I've found a working one:

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

where lcm is least common multiple. Does anyone have a better one?

Note

I'm using python, I think code should be readable even to who doesn't know python, if you find anything that is not clear just ask in the comments. I'm interested in the algorithm and the math, the idea.

解决方案

[Edit] One last thought, which (IMO) is important enough that I'll put it at the beginning: if you're collecting a bunch of totients at once, you can avoid a lot of redundant work. Don't bother starting from large numbers to find their smaller factors -- instead, iterate over the smaller factors and accumulate results for the larger numbers.

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

This takes just 8ms on my desktop.

---

The Wikipedia page on the Euler totient function has some nice mathematical results.

counts the numbers coprime to and smaller than each divisor of : this has a trivial* mapping to counting the integers from to , so the sum total is .

* by the second definition of trivial

This is perfect for an application of the Möbius inversion formula, a clever trick for inverting sums of this exact form.

This leads naturally to the code

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

There exist better implementations of the Möbius function, and it could be memoized for speed, but this should be easy enough to follow.

The more obvious computation of the totient function is

In other words, fully factor the number into unique primes and exponents, and do a simple multiplication from there.

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

Again, there exist better implementations of prime_factors, but this is meant for easy reading.

---

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

Using this code to calculate the totients of all numbers from 1 to 9999 on my desktop, averaging over 5 runs,

• totient1 takes forever
• totient2 takes 10s
• totient3 takes 1.3s
• totient4 takes 1.3s

这篇关于N 以下有多少个数是 N 的互质数?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！