Published: 2018-09-04

最大公约数算法

1 问题

两个不全为0的非负整数m和n的最大公约数记做 gcd(m,n),代表能够整除(即余数为0)m和n的最大正整数。 求gcd(m,n)的值。

2 欧几里得算法

思路: gcd(m,n) = gcd(n, m mod n ) (m mod n 表示m除以n之后的余数) 。因为gcd(m,0) = m, 所以m最后的取值就是m和n初值的最大公约数

def gcd(m,n):

    while n != 0:
        r = m % n
        m = n
        n = r
    return m

3 连续整数检测算法

思路: gcd(m,n)的必然不能大于min(m,n), 从t={m,n}开始,检查t能不能整除m和n,如果能则t是最大公约数,如果不能则t=t-1重复这个步骤, 直到t不大于0。

def gcd(m,n):

    # corner case
    if m == 0: return n
    if n == 0: return m

    t = min(m, n)
    while t > 0 and not (m % t == 0 and n % t == 0):
        t -= 1
    return t

4 中学时计算最大公约数

思路: 找到m和n的所有质因数,找出两者质因数公共的部分,这些质因数相乘即是最大公约数。

4.1 埃拉托色尼筛选法

衍生问题:如何产生一个不大于给定正整数n的连续质数序列。

思路:一开始初始化一个2-n的连续整数序列作为候选质数,然后在第一个循环中消除2的倍数, 第二个循环中消除3的倍数,第三个循环中消除5的倍数(4已结消除了)。 这样一直下去直到序列中没有可消除的元素为止。

一个优化:考虑到正在消除p的倍数,第一个值得考虑的倍数是p*p, 因为2p,3p,…,(p-1)p已经在先前的步骤中消去。所以如果p*p>n,那么自然没必要继续消除下去

4.2 实现

import math
from functools import reduce


def find_primes(n):
    "埃拉托色尼筛选法, 找出不超过n的质数"

    l = list(range(2, n+1))

    # 心得: 这里遇到在迭代列表的时候,要移除列表的元素,但移除完后列表的后续索引就变了,索引采用了一个额外的数组来表示删除与否
    # l = [1,2,3,4] l_index=[True, True, False, True], 就代表着 [1,2,4]
    l_index = [True] * (n-1)

    p = math.floor(math.sqrt(n))

    for i in range(0, p):
        if not l_index[i]:
            continue
        for j in range(i+1, n-1):
            if not l_index[j]:
                continue
            if l[j] % l[i] == 0:
                l_index[j] = False

    return [ l[i] for i in range(0, n-1) if l_index[i]]


def find_prime_factors(n):
    "找出n的质因数"

    result = []
    all_primes = find_primes(n)
    while n >1:
        for p in all_primes:
            if n % p == 0:
                result.append(p)
                n = int(n / p)
    return result


def get_common_part(l1, l2):
    "获取两个数字列表共同的元素"

    # 先排序再比较
    l1 = sorted(l1)
    l2 = sorted(l2)
    l1_len = len(l1)
    l2_len = len(l2)
    i, j = 0, 0
    result = []

    while i < l1_len and j < l2_len:
        if l1[i] == l2[j]:
            result.append(l1[i])
            i += 1
            j += 1
        elif l1[i] > l2[j]:
            j += 1
        else:
            i += 1
    return result




def gcd(m, n):
    if m == 0: return n
    if n == 0: return m

    common_part = get_common_part(*map(find_prime_factors, [m, n]))

    if not common_part:
        return 1

    return reduce(lambda x, y: x*y, common_part)

Author: Nisen

Email: imnisen@163.com

Emacs 25.2.1 (Org mode 8.2.10)