找出小于N的最大平方根约数

Question

实际上，给定一个可能非常大的偶数N，我想找到N = F * R，其中gcd（F，R）= 1，F> R，并且F尽可能小（因为我将完全分解F）。问题的核心是找到最大的除数R，其中R < sqrt（N）。

例如，N = 36应该给出F = 9和R = 4。请注意，R不一定是质数或质数幂。请注意，我没有将N分解。对F和R的唯一限制是它们是互质的。

这是我的快速和天真版本，它正在工作：

def factor_partial(N):
    for R in xrange(int(math.sqrt(N)),1,-1):
        if N%R == 0 and gcd(R,N/R) == 1:
            return N/R, R

我想到的另一种方法是按照递增顺序找到除数，并在此过程中删除任何非除数的倍数。类似于：

def factor_partial(N):
    i = range(2, int(sqrt(N)) + 1)
    while i:
        if N % i[0] != 0:
            remove_multiples(i, i[0]) #without removing i[0]
        else:
            if gcd(i[0], N/i[0]) == 1:
                R = i[0]
        i.pop(0) #remove i[0]

    return N/R, R

我认为这将会很慢且内存占用大，但如果i是一个生成器，那么它可能会更高效。我还没有多少使用生成器的经验。

我能改进第一个版本吗？第二个版本可行吗（我该如何做）？是否有完全不同的方法更好？

寻找Python、C或伪代码的答案。

---

这是一个关于数论课程的项目。我正在实现一个基于Pocklington的素性测试。虽然我需要一个因数分解算法，但我们还没有学习过任何算法，我可能不会使用像二次筛法这样超出课程范围的算法。我正在寻求针对所提出问题的具体帮助。

- jmilloy

---

你的第二种方法是Eratosthenes筛法。我建议你也看一下这个问题。 - Blender

@Blender 我为什么需要存储前面的数字？请注意，我不是在寻找质因数，实际上这段代码是筛法的一种补充，因为我正在删除i的倍数，而i不是N的除数。我只关心最大的（质数或复合数）除数，我将其保存为R，因此在检查完每个数字后都会丢弃它。 - jmilloy

让我们在聊天中继续这个讨论：http://chat.stackoverflow.com/rooms/5352/discussion-between-jmilloy-and-blender - jmilloy

这是一个使用代码解释的数学/算法问题，对我来说更像是一个数学问题，而不是编程问题。 - razlebe

@razlebe 是否有更好的StackExchange网站可以提问这种算法问题？ - jmilloy

Answer 1

维基百科有一个很好的因数分解算法列表: http://en.wikipedia.org/wiki/Integer_factorization#Factoring_algorithms

你的第二种方法有效地使用了筛子，并且当N是一些小质数的倍数时，具有快速减少问题的好特性。通过循环素数而不是所有可能的2..sqrt(n)除数来改进代码。

此外，您可能希望首先进行素性测试，以便在执行其他工作之前知道N是合成数。

您的笔记说您没有因式分解N，但问题是相关的。寻找F和R相当于探索N的质因数的非重叠组合。

对于N==36的情况，N的质因数分解为2、2、3、3。F和R的因子必须包括所有这些（以便F*R==N），并且不能重叠（以便GCD(F,R)==1）。 因此4和9立即出现。

一个更具指导性的例子可能是N==23256。它的因式分解是2,2,2,3,3,17,19。由于F和R之间不能有重叠，每个质数基数只能进入两个桶中的一个（即你要么得到所有的二，要么没有）。因此，我们可以将这些因子分组为8,9,17,19。为了找到R，我们希望选择那些因子的组合尽可能大但不超过152.49，即23256的平方根。我们的选择是{8}，{9}，{8,9}，{8,17}，{8,19}。其中最大的是8*19，等于152。相应的F是17*19或153。

上述choices是通过[choice for choice in powerset([8,9,17,19]) if prod(choice) < math.sqrt(N)]计算出来的。

因此，整个程序基本上归结为：

prime_factors = factorize(N)      # [2,2,2,3,3,17,19]
clusters = [p**e for p, e in collections.Counter(prime_factors).items()]  # [8,9,17,19]
R = max(prod(group) for group in powerset(clusters) if prod(group) < math.sqrt(N))
F = N // R

当powerset搜索超过N的平方根时，可以通过修剪集合的生成来加快搜索速度。

请记住，分解是计算密集型的，而幂集增长非常快，但它很可能比从N的平方根开始进行许多除法的原始算法要少得多。

- Raymond Hettinger

---

假设N为36，则我要寻找的解为R=4，F=9。如果我循环质数，我找不到R=4。如果我们限制只检查质数，如何恢复此解？ - jmilloy

“4”在因数分解中根本不应被列出，你需要找到的是质因数。在这种情况下，你应该找到“2”两次。 - jsbueno

@jsbueno 不是的。我没有分解N。F和R唯一的限制是它们互质。 - jmilloy

@jmilloy：重点是解决问题的最简单方法是先分解N，然后将因子分组。 - Chris Dodd

谢谢。这有点滑稽 - 我正在编写一个素性测试，重点是它应该比分解N更容易找到N-1 = F * R并分解F。我决定尝试最小化F，结果发现（似乎）这与一开始分解N-1一样困难。 - jmilloy

Answer 2

你能否对N进行质因数分解，然后找到所有质因子的最大乘积组合，使其小于sqrt(N)？

例如，对于36，它会发现质因数分解为2*2*3*3。然后您将尝试所有不同的质数组合：

2 = 2
3 = 3
2*2 = 4
2*3 = 6
3*3 = 9
2*2*3 = 12
2*3*3 = 18
2*2*3*3 = 36

您知道这些都是36的因数，因此您找到最大的一个，使其小于sqrt(36)，结果是4。

但是，除非您已经有一个现有的质数列表或质因数分解，或者拥有一些令人惊叹的质因数分解算法，或者正在处理极大的数字，否则我真的看不出来它如何比您的第一个版本更快。

但即使如此（回到第一个版本的狂热），O(sqrt(n))也是一个相当快的运行时间，只需要O(1)的内存，所以第一个算法可能只是最好的选择。我不认为它会很慢，特别是在现代计算机上使用C语言。

- Andrew Rasmussen

---

嗯，是的，N的分解能解决整个问题...但这个算法是代替获取N的分解的 :P - jmilloy

哦哈哈。那么我会坚持你的第一个版本。它在内存和运行时间方面都非常高效。祝你好运！ - Andrew Rasmussen

不需要考虑像 R=23 这样的组合，因为 F 就会得到另外的质数 F=23，这些质数就不是相对质数了。所以在给定的因子中，你要么得到所有的二，要么一个也没有。因此，可能的候选因子是 {1, 4, 9, 36}。 - Raymond Hettinger

很抱歉，我不太明白你的解释？ - Andrew Rasmussen

OP想要找到两个因数F和R，它们的乘积为N，其中R尽可能大（最大到sqrt(N)），并且它们互质。后面这个限制意味着2*3不是R的候选项（因为相应的F=N/R不与R互质）。这个限制有助于大大减少搜索空间。 - Raymond Hettinger

Answer 3

def factor_partial(N):
    R = int(sqrt(float(N)))
    sieve = [1, 1] + [0] * (R-1)
    for i in xrange(2, R) :
        if sieve[i]==0 :
            j=i*i;
            while j<=R :
                sieve[j]=1
                j = j + i
    primes = [i for i in xrange(R+1) if sieve[i]==0]

    saveN = N
    primepower_divisors = []
    for i in primes :
        if N%i == 0 :
            term = 1
            while N%i == 0 :
                term = term * i
                N = N / i
            primepower_divisors = primepower_divisors + [term]
            if N==1 : break

    largecoprime_divisors = [1]
    for i in primepower_divisors :
        for j in largecoprime_divisors :
            largecoprime_divisors = largecoprime_divisors + [i*j]

    F = min([i for i in largecoprime_divisors if i>R])
    return F, saveN/F

我使用筛法计算质数列表（在计算质数列表时有很多优化方法可用）。 我们可以利用这个事实，即不存在质数p使得F%p == 0且R%p == 0。 因为gcd(F, R)=1。

- anshul410

---

通常情况下，min([i for i in largecoprime_divisors if i>R]) 会失败，因为它的序列为空。有时候它比我的第一个版本更慢。 - jmilloy

是的，这个速度比较慢.. 我刚刚发现了.. :) 也许你应该使用Polard的Rho.. - anshul410