找到整数的所有因子的有效算法是什么？

Question

我正在编写一个非常简单的程序来检查一个数字是否可以整除另一个数字：

// use the divider squared to reduce iterations
for(divider = 2; (divider * divider) <= number; divider++)
    if(number % divider == 0)
        print("%d can divided by %d\n", number, divider);

现在我很好奇这项任务是否可以通过找到数字的平方根并将其与除法器进行比较来完成。但是，似乎 sqrt() 并不能真正提高效率。如何在 C 中处理 sqrt() 以及如何提高 sqrt() 的效率？另外，有没有其他方法可以更有效地解决问题？

还有，

number % divider == 0

用于测试除法器是否可以均匀地除数，除了使用％之外还有没有更有效的测试方法？

Answer 1

我不会讨论找到整数所有因子的最佳算法是什么。相反，我想对您当前的方法发表评论。

有条件测试用例需要考虑

1. (divider * divider) <= number
2. divider <= number/divider
3. divider <= sqrt(number)

查看Conditional tests in primality by trial division了解更多详情。

使用的案例取决于您的目标和硬件。

案例 1 的优点是不需要除法。但是，当divider*divider 大于最大整数时，它可能会溢出。情况二没有溢出问题，但需要除法。对于 case3，sqrt 只需要计算一次，但它要求 sqrt 函数得到正确的完美平方。

但是还有一些其他的事情需要考虑很多指令集，包括 x86 指令集，在进行除法时也会返回余数。由于您已经在使用number % divider，这意味着您在使用number / divider 时可以免费获得它。

因此，情况 1 仅适用于除法和余数不是在一条指令中计算且您不必担心溢出的系统。

在案例 2 和案例 3 之间，我认为主要问题还是指令集。如果 sqrt 与 case2 相比太慢，或者您的 sqrt 函数无法正确计算完全平方，请选择 case 2。如果指令集没有在一条指令中计算除数和余数，则选择案例 3。

对于 x86 指令集案例 1，案例 2 和案例 3 应该提供基本相同的性能。所以 should 没有理由使用案例 1（但是请参阅下面的一个微妙点）。 C 标准库保证完美正方形的sqrt 正确完成。所以情况3也没有缺点。

但是关于案例 2 有一个微妙的点。我发​​现一些编译器不承认除法和余数是一起计算的。例如在下面的代码中

for(divider = 2; divider <= number/divider; divider++)
    if(number % divider == 0)

GCC 生成两个除法指令，即使只需要一个。解决此问题的一种方法是像这样保持划分和提醒关闭

divider = 2, q = number/divider, r = number%divider
for(; divider <= q; divider++, q = number/divider, r = number%divider)
    if(r == 0)

在这种情况下，GCC 只产生一个除法指令，并且 case1、case 2 和 case 3 具有相同的性能。但是这段代码的可读性比

int cut = sqrt(number);
for(divider = 2; divider <= cut; divider++)
    if(number % divider == 0)

所以我认为总体情况 3 是最好的选择，至少对于 x86 指令集而言。

Answer 2

但是，sqrt() 似乎并不能真正提高效率

这是意料之中的，因为每次迭代节省的乘法主要由循环内慢得多的除法运算支配。

另外，number % divider = 0是用来测试除法器是否能均匀除数，除了使用%之外，还有没有更有效的测试方法？

我不知道。检查a % b == 0 是否至少与检查a % b = c 的某些c 一样困难，因为我们可以使用前者来计算后者（加上一个额外的加法）。至少在英特尔架构上，计算后者的计算成本与除法一样昂贵，这是典型现代处理器中最慢的运算之一。

如果您希望显着获得更好的性能，则需要更好的分解算法，其中there are plenty。一个特别简单的运行时 O(n^1/4) 是Pollard's ρ algorithm。你可以找到一个简单的 C++ 实现in my algorithms library。对 C 的适应留给读者作为练习：

int rho(int n) { // will find a factor < n, but not necessarily prime
  if (~n & 1) return 2;
  int c = rand() % n, x = rand() % n, y = x, d = 1;
  while (d == 1) {
    x = (1ll*x*x % n + c) % n;
    y = (1ll*y*y % n + c) % n;
    y = (1ll*y*y % n + c) % n;
    d = __gcd(abs(x - y), n);
  }
  return d == n ? rho(n) : d;
}

void factor(int n, map<int, int>& facts) {
  if (n == 1) return;
  if (rabin(n)) { // simple randomized prime test (e.g. Miller–Rabin)
    // we found a prime factor
    facts[n]++;
    return;
  }
  int f = rho(n);
  factor(n/f, facts);
  factor(f, facts);
}

从其素因子构造n 的因子是一件容易的事。只需将所有可能的指数用于找到的质因数并以每种可能的方式组合它们。

Answer 3

在 C 中，您可以使用标头 <math.h> 中的 sqrt() 系列函数对浮点数求平方根。

取平方根通常比除法要慢，因为取平方根的算法比除法算法更复杂。这不是 C 语言的属性，而是执行程序的硬件的属性。在现代处理器上，求平方根可以和除法一样快。例如，这适用于 Haswell 微架构。

但是，如果算法改进良好，sqrt() 调用的速度稍慢通常无关紧要。

要只比较 number 的平方根，请使用如下代码：

#include <math.h>

/* ... */

int root = (int)sqrt((double)number);
for(divider = 2; divider <= root; divider++)
    if(number % divider = 0)
        print("%d can divided by %d\n", number, divider);

Answer 4

这只是我的一个想法，如果有错误请评论和批评。

这个想法是预先计算一定范围内的所有素数并将其用作表格。

循环遍历表，检查素数是否是一个因子，如果是，则增加该素数的计数器，如果不是，则增加索引。当索引到达末尾或要检查的素数超过输入时终止。

最后，结果是输入的所有主要因子及其计数的表格。那么生成所有自然因素应该是微不足道的，不是吗？

最坏的情况，循环需要走到尽头，然后需要6542次迭代。

考虑到输入是 [0, 4294967296] 这类似于O(n^3/8)。

这是实现此方法的 MATLAB 代码：

如果p 由p=primes(65536); 生成，则此方法适用于[0, 4294967296] 之间的所有输入（但未经测试）。

function [ output_non_zero ] = fact2(input, p)
    output_table=zeros(size(p));
    i=1;
    while(i<length(p));
        if(input<1.5)
            break; 
            % break condition: input is divided to 1,
            % all prime factors are found. 
        end
        if(rem(input,p(i))<1)
            % if dividable, increament counter and don't increament index
            % keep checking until not dividable
            output_table(i)=output_table(i)+1;
            input = input/p(i);
        else
            % not dividable, try next
            i=i+1;
        end
    end

    % remove all zeros, should be handled more efficiently
     output_non_zero = [p(output_table~=0);...
                        output_table(output_table~=0)];
     if(input > 1.5)
         % the last and largest prime factor could be larger than 65536
         % hence would skip from the table, add it to the end of output
         % if exists
         output_non_zero = [output_non_zero,[input;1]];
     end
end

测试

p=primes(65536);
t = floor(rand()*4294967296);
b = fact2(t, p);
% check if all prime factors adds up and they are all primes
assert((prod(b(1,:).^b(2,:))==t)&&all(isprime(b(1,:))), 'test failed');