3

我知道已经有一个与此类似的问题，但我想使用 GMPY2（或与 GMP 类似的东西）加速它。这是我当前的代码，它很不错，但它可以更好吗？

编辑：新代码，检查除数 2 和 3

def factors(n):
    result = set()
    result |= {mpz(1), mpz(n)}


    def all_multiples(result, n, factor):
        z = mpz(n)
        while gmpy2.f_mod(mpz(z), factor) == 0:
            z = gmpy2.divexact(z, factor)
            result |= {mpz(factor), z}
        return result

    result = all_multiples(result, n, 2)
    result = all_multiples(result, n, 3)

    for i in range(1, gmpy2.isqrt(n) + 1, 6):
        i1 = mpz(i) + 1
        i2 = mpz(i) + 5
        div1, mod1 = gmpy2.f_divmod(n, i1)
        div2, mod2 = gmpy2.f_divmod(n, i2)
        if mod1 == 0:
            result |= {i1, div1}
        if mod2 == 0:
            result |= {i2, div2}
    return result

如果可能的话，我也对仅包含除数的实现感兴趣n^(1/3) and 2^(2/3)*n(1/3)

例如，mathematicafactor()比 python 代码快得多。我想分解 20 到 50 个十进制数字之间的数字。我知道 ggnfs 可以在不到 5 秒的时间内解决这些问题。

如果python中也存在任何实现快速分解的模块，我很感兴趣。

4

2 回答 2

7

我刚刚对您的代码进行了一些快速更改，以消除多余的名称查找。该算法仍然相同，但在我的计算机上它的速度大约是原来的两倍。

import gmpy2
from gmpy2 import mpz

def factors(n):
    result = set()
    n = mpz(n)
    for i in range(1, gmpy2.isqrt(n) + 1):
        div, mod = divmod(n, i)
        if not mod:
            result |= {mpz(i), div}
    return result

print(factors(12345678901234567))

其他建议将需要有关数字大小等的更多信息。例如，如果您需要所有可能的因素，那么从所有主要因素中构建这些可能会更快。这种方法可以让您在继续进行时减少范围语句的限制，并且还可以让您增加 2（在删除 2 的所有因素之后）。

更新 1

我对您的代码进行了一些额外的更改。我不认为你的all_multiplies()功能是正确的。您的range()陈述不是最佳的，因为 2 再次检查，但我的第一次修复使它变得更糟。

新代码延迟计算余数，直到它知道余数为 0。我还尝试尽可能多地使用内置函数。例如，mpz % integer比 gmpy2.f_mod(mpz, integer) 或 gmpy2.f_mod(integer, mpz) 更快，其中integer是一个普通的 Python 整数。

import gmpy2
from gmpy2 import mpz, isqrt

def factors(n):
    n = mpz(n)

    result = set()
    result |= {mpz(1), n}

    def all_multiples(result, n, factor):
        z = n
        f = mpz(factor)
        while z % f == 0:
            result |= {f, z // f}
            f += factor
        return result

    result = all_multiples(result, n, 2)
    result = all_multiples(result, n, 3)

    for i in range(1, isqrt(n) + 1, 6):
        i1 = i + 1
        i2 = i + 5
        if not n % i1:
            result |= {mpz(i1), n // i1}
        if not n % i2:
            result |= {mpz(i2), n // i2}
    return result

print(factors(12345678901234567))

我会改变你的程序，只找到所有小于平方根的素因子n，然后再构造所有的协因子。然后n每次找到一个因子时都会减少，检查是否n是素数，如果不是素数，则只寻找更多的因子n。

更新 2

pyecm模块应该能够分解您尝试分解的大小数字。以下示例在大约一秒钟内完成。

>>> import pyecm
>>> list(pyecm.factors(12345678901234567890123456789012345678901, False, True, 10, 1))
[mpz(29), mpz(43), mpz(43), mpz(55202177), mpz(2928109491677), mpz(1424415039563189)]
于 2014-05-17T13:50:16.973 回答
0

Internet 中存在不同的 Python 分解模块。但是，如果您想自己实现因式分解（不使用外部库），那么我可以建议您非常快速且非常容易地实现Pollard-Rho Algorithm。我在下面的代码中完全实现了它，如果您不想阅读，您只需直接向下滚动到我的代码（在答案的底部）。

Pollard-Rho 算法很有可能在 的时间内找到最小的非平凡因子P（不等于1或N）O(Sqrt(P))。为了比较，您在问题中实现的Trial DivisionO(P)算法需要时间才能找到 factor P。这意味着例如，如果一个素数因子P = 1 000 003，则试除法将在1 000 003除法运算之后找到它，而 Pollard-Rho 平均会在1 000运算 ( Sqrt(1 000 003) = 1 000) 之后找到它，这要快得多。

为了使 Pollard-Rho 算法更快，我们应该能够检测素数，将它们排除在因式分解之外并且不要等待不必要的时间，因为在我的代码中我使用了Fermat Primality Test，它非常快速且易于实现7-9行代码。

Pollard-Rho 算法本身很短，13-15 行代码，你可以在我的函数的最底部看到它pollard_rho_factor()，其余代码行是补充的辅助函数。

我从头开始实现了所有算法，而不使用额外的库（random模块除外）。这就是为什么你可以在那里看到我的gcd()函数，尽管你可以使用内置的 Python 的math.gcd()代替（它找到最大公约数）。

您可以在我的代码中看到函数Int()，它仅用于将 Python 的整数转换为GMPY2。GMPY2 ints 将使算法更快，您可以使用 Pythonint(x)代替。我没有使用任何特定的 GMPY2 函数，只是将所有整数转换为 GMPY2 整数，以实现大约 50% 的加速。

例如，我考虑了Pi的前 190 位数字！将它们分解需要 3-15 秒。Pollard-Rho 算法是随机的，因此每次运行需要不同的时间来分解相同的数字。您可以再次重新启动程序，看看它会打印不同的运行时间。

当然，分解时间很大程度上取决于素数除数的大小。一些 50-200 位数字可以在几分之一秒内分解，有些需要几个月。我的示例 190 位 Pi 的素因数非常小，除了最大的一个，这就是它速度快的原因。Pi 的其他数字可能没有那么快分解。所以数字的位数并不重要，只有素数的大小很重要。

我有意pollard_rho_factor()将函数实现为一个独立的函数，而不是将其分解为更小的独立函数。虽然它打破了 Python 的风格指南，但（我记得）它建议不要嵌套函数并将所有可能的函数放在全局范围内。样式指南还建议在脚本的第一行中在全局范围内进行所有导入。我特意做了一个函数，以便它可以轻松复制粘贴并完全准备好在您的代码中使用。费马素数测试is_fermat_probable_prime()子功能也是可复制粘贴的，无需额外依赖即可工作。

在极少数情况下，Pollard-Rho 算法可能无法找到非平凡素因子，尤其是对于非常小的因子，例如，您可以用少量替换ninside ，然后发现 Pollard-Rho 失败。对于如此小的失败因素，您可以轻松地使用您在问题中实现的Trial Division算法。test()4

在线尝试！

def pollard_rho_factor(N, *, trials = 16):
    # https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    import math, random
    
    def Int(x):
        import gmpy2
        return gmpy2.mpz(x) # int(x)
    
    def is_fermat_probable_prime(n, *, trials = 32):
        # https://en.wikipedia.org/wiki/Fermat_primality_test
        import random
        if n <= 16:
            return n in (2, 3, 5, 7, 11, 13)
        for i in range(trials):
            if pow(random.randint(2, n - 2), n - 1, n) != 1:
                return False
        return True
    
    def gcd(a, b):
        # https://en.wikipedia.org/wiki/Greatest_common_divisor
        # https://en.wikipedia.org/wiki/Euclidean_algorithm
        while b != 0:
            a, b = b, a % b
        return a
        
    def found(f, prime):
        print(f'Found {("composite", "prime")[prime]} factor, {math.log2(f):>7.03f} bits... {("Pollard-Rho failed to fully factor it!", "")[prime]}')
        return f
        
    N = Int(N)
    
    if N <= 1:
        return []
    
    if is_fermat_probable_prime(N):
        return [found(N, True)]
    
    for j in range(trials):
        i, stage, y, x = 0, 2, Int(1), Int(random.randint(1, N - 2))
        while True:
            r = gcd(N, abs(x - y))
            if r != 1:
                break
            if i == stage:
                y = x
                stage <<= 1
            x = (x * x + 1) % N
            i += 1
        if r != N:
            return sorted(pollard_rho_factor(r) + pollard_rho_factor(N // r))
    
    return [found(N, False)] # Pollard-Rho failed

def test():
    import time
    # http://www.math.com/tables/constants/pi.htm
    # pi = 3.
    #     1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679
    #     8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196
    # n = first 190 fractional digits of Pi
    n =   1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489
    tb = time.time()
    print('N:', n)
    print('Factors:', pollard_rho_factor(n))
    print(f'Time: {time.time() - tb:.03f} sec')
    
test()

输出：

N: 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489
Found prime factor,   1.585 bits...
Found prime factor,   6.150 bits...
Found prime factor,  20.020 bits...
Found prime factor,  27.193 bits...
Found prime factor,  28.311 bits...
Found prime factor, 545.087 bits...
Factors: [mpz(3), mpz(71), mpz(1063541), mpz(153422959), mpz(332958319), mpz(122356390229851897378935483485536580757336676443481705501726535578690975860555141829117483263572548187951860901335596150415443615382488933330968669408906073630300473)]
Time: 2.963 sec
于 2021-10-21T04:59:20.037 回答