【问题标题】:Most efficient way to find all factors with GMPY2 (or GMP)?使用 GMPY2(或 GMP)查找所有因素的最有效方法?
【发布时间】:2014-05-17 05:46:32
【问题描述】:

我知道已经有一个与此类似的问题,但我想使用 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) 内使用除数的实现感兴趣

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

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

【问题讨论】:

  • 答案很大程度上取决于n的大小。你能告诉我们你想考虑多大的 n(多少个十进制数字)?

标签: python factorization gmpy


【解决方案1】:

我刚刚对您的代码进行了一些快速更改,以消除多余的名称查找。算法还是一样的,但在我的电脑上速度是原来的两倍。

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

【讨论】:

    【解决方案2】:

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

    Pollard-Rho 算法很有可能在O(Sqrt(P)) 的时间内找到最小的非平凡因子P(不等于1N)。为了比较,您在问题中实现的Trial Division 算法需要O(P) 时间才能找到因子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() 代替(它找到 Greatest Common Divisor)。

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

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

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

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

    在极少数情况下,Pollard-Rho 算法可能无法找到非平凡素因子,尤其是对于非常小的因子,例如,您可以将 test() 内的 n 替换为小数字 4 并看到 Pollard-罗失败了。对于如此小的失败因素,您可以轻松使用您在问题中实现的 Trial Division 算法。

    Try it online!

    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
    

    【讨论】:

    • 这是我的一个非常古老的问题。之后,我确实尝试用 C++ 等更快的语言实现 Pollard rho 和二次筛等算法。但我很欣赏这个答案。
    • @qwr 我还在 Python 和 C++ 中实现了二次筛。二次筛真的非常快,即使它只有两个大小相等的因子,它也可以在 10 分钟内分解 256 位。当然还有更快的 GNFS 算法,但实现起来非常复杂。二次筛具有中等的速度和复杂性。在这里,我实现了 Pollard Rho,作为一个实现起来非常简单但速度很快的算法的例子。 Quadratic Sieve 是一种在 StackOverflow 上发布的复杂方法。 StackOverflow 的帖子有 50 000 个字符的限制,并且二次筛的来源更大。
    • @qwr 我非常喜欢二次筛。你能说出你的 C++ 二次筛可以分解的最大位数是多少吗?它花费了多少时间来做到这一点?我想比较我的时间,以确定我是否有优化版本或者它需要更多优化。如果您不记得时间或没有时间进行实验,请跳过此评论请求。
    • 我很久以前写的代码从来都不是很漂亮,那里有很多完全浪费的东西,但是你可以看到它github.com/jxu/Factorize
    猜你喜欢
    • 2020-10-02
    • 1970-01-01
    • 2011-10-11
    • 2010-09-29
    • 2021-11-11
    • 1970-01-01
    • 2015-03-25
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多