整数分解

Question

在回答另一个问题时，我偶然发现了如何不使用Symbolic Math Toolbox找到整数的所有因子。

例如：

factor(60)

返回：

 2     2     3     5

---

unique(factor(60))

因此会返回所有质数因子，"1" 除外。

 2     3     5

---

我正在寻找一个函数，它将返回所有因数（1和该数字本身不重要，但它们很好）

当x = 60时，预期输出结果为：

 1     2     3     4     5     6    10    12    15    20    30    60     

我想出了这个相当笨重的解决方案，除此之外它可能可以进行向量化处理，难道没有更优雅的解决方案吗？

---

x = 60;

P = perms(factor(x));
[n,m] = size(P);
Q = zeros(n,m);
for ii = 1:n
    for jj = 1:m
        Q(ii,jj) = prod(P(ii,1:jj));
    end
end

factors = unique(Q(:))'

我认为，这个解决方案对于某些大数字将会失败，因为 perms 要求向量长度小于11。

- Robert Seifert

---

你的代码没有生成12和20，仅供参考。 - Mahm00d

@Mahm00d：你说得对，甚至更糟。我已经纠正了它。 - Robert Seifert

Answer 1

通过将一个包含整数1到n的向量除以一个数字n，然后找到除以1余数为零（即整数结果）的位置，就可以找到数字n的所有因子：

>> n = 60;
>> find(rem(n./(1:n), 1) == 0)

ans =

     1     2     3     4     5     6    10    12    15    20    30    60

- gnovice

---

我之前考虑过类似的东西，但从未听说过“rem”。这非常优雅！ - Robert Seifert

不需要除法，仅使用“rem”即可（请参见我的答案）。 - Luis Mendo

这种暴力方法对于较大的整数来说速度很慢，更不用说它很容易抛出内存错误（例如n=1e9需要约8GB的内存）。 - Amro

@Amro：对于更大的整数，内存确实是一个问题，但可以通过一些技巧来缓解（其中一些在您的答案中提到）。我只是给出了一个快速而简单的答案，适用于大多数情况，并且能够帮助OP。 :) - gnovice

很难决定哪个答案是“最好的”，对于小数字，你/路易斯的方法非常有建设性，但对于相当大的数字，阿姆罗的方法是无与伦比的，这是非常深思熟虑的。 - Robert Seifert

@thewaywewalk 完全可以理解。别担心。 :) - gnovice

Answer 2

以下是寻找整数因子的六种不同实现方法的比较:

function [t,v] = testFactors()
    % integer to factor
    %{45, 60, 2059, 3135, 223092870, 3491888400};
    n = 2*2*2*2*3*3*3*5*5*7*11*13*17*19;

    % functions to compare
    fcns = {
        @() factors1(n);
        @() factors2(n);
        @() factors3(n);
        @() factors4(n);
        %@() factors5(n);
        @() factors6(n);
    };

    % timeit
    t = cellfun(@timeit, fcns);

    % check results
    v = cellfun(@feval, fcns, 'UniformOutput',false);
    assert(isequal(v{:}));
end

function f = factors1(n)
    % vectorized implementation of factors2()
    f = find(rem(n, 1:floor(sqrt(n))) == 0);
    f = unique([1, n, f, fix(n./f)]);
end

function f = factors2(n)
    % factors come in pairs, the smaller of which is no bigger than sqrt(n)
    f = [1, n];
    for k=2:floor(sqrt(n))
        if rem(n,k) == 0
            f(end+1) = k;
            f(end+1) = fix(n/k);
        end
    end
    f = unique(f);
end

function f = factors3(n)
    % Get prime factors, and compute products of all possible subsets of size>1
    pf = factor(n);
    f = arrayfun(@(k) prod(nchoosek(pf,k),2), 2:numel(pf), ...
        'UniformOutput',false);
    f = unique([1; pf(:); vertcat(f{:})])'; %'
end

function f = factors4(n)
    % http://rosettacode.org/wiki/Factors_of_an_integer#MATLAB_.2F_Octave
    pf = factor(n);                    % prime decomposition
    K = dec2bin(0:2^length(pf)-1)-'0'; % all possible permutations
    f = ones(1,2^length(pf));
    for k=1:size(K)
      f(k) = prod(pf(~K(k,:)));        % compute products 
    end; 
    f = unique(f);                     % eliminate duplicates
end

function f = factors5(n)
    % @LuisMendo: brute-force implementation
    f = find(rem(n, 1:n) == 0);
end

function f = factors6(n)
    % Symbolic Math Toolbox
    f = double(evalin(symengine, sprintf('numlib::divisors(%d)',n)));
end

结果如下：

>> [t,v] = testFactors();
>> t
t =
    0.0019        % factors1()
    0.0055        % factors2()
    0.0102        % factors3()
    0.0756        % factors4()
    0.1314        % factors6()

>> numel(v{1})
ans =
        1920

尽管第一个矢量化版本是最快的，但由于自动JIT优化，等效的基于循环的实现（factors2）也不落后太多。

请注意，我必须禁用暴力实现（factors5()），因为它会抛出内存不足错误（将1:3491888400存储为双精度需要超过26GB的内存！）。 对于大整数，无论从空间还是时间上来看，这种方法显然都不可行。

结论：请使用以下矢量化实现 :)

n = 3491888400;
f = find(rem(n, 1:floor(sqrt(n))) == 0);
f = unique([1, n, f, fix(n./f)]);

- Amro

---

不错的想法！你知道使用mod而不是rem可以将时间减半吗？（这就是为什么我从未听说过rem，我总是使用mod，在大多数情况下它们实际上是相同的）- PS：我更愿意把你的结论放在最前面，我不知道有多少人关心其他5种方法;) - Robert Seifert

@thewaywewalk：实际上，我并没有看到rem和mod之间的时间差异，正如你所说，它们在列出的特殊情况以外几乎是相同的。我正在运行最新的64位Windows MATLAB版本。 - Amro

我用你的计时函数再次测试了一下，结果要么两者时间相同，要么（在大多数情况下）其中一个快了5倍。不过，这并不重要。 - Robert Seifert

Answer 3

@gnovice的回答有所改进，可以避免使用除法操作：仅用rem即可：

n = 60;
find(rem(n, 1:n)==0)

- Luis Mendo

---

好的发现！我完全没注意到。 - gnovice

我认为保持gnovices的答案被选为接受是公平的，因为他有这个想法;) 但无论如何+1！ - Robert Seifert