Факторизация целого числа

Question

Отвечая на другое, я наткнулся на вопрос, каким образом я мог бы найти все факторы целого числа без набора инструментов Symbolic Math.

Например:

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.

Answer 1

Вот сравнение шести различных реализаций для нахождения факторов целого числа:

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

Хотя первая векторизованная версия является самой быстрой, эквивалентная реализация на основе цикла (factors2) не сильно отстает, благодаря автоматической оптимизации JIT.

Обратите внимание, что мне пришлось отключить реализацию грубой силы (factors5()) потому что он выдает ошибку нехватки памяти (сохраняя вектор 1:3491888400 в двойной точности требуется более 26 ГБ памяти!). Этот метод, очевидно, не подходит для больших целых чисел ни в пространстве, ни во времени.

Вывод: используйте следующую векторизованную реализацию:)

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

Answer 2

Вы можете найти все факторы ряда n разделив его на вектор, содержащий целые числа от 1 до nзатем найти остаток от деления на 1, равный нулю (т. е. целочисленные результаты):

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

ans =

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

Answer 3

Улучшение по сравнению с ответом @gnovice - пропустить операцию деления: rem одного достаточно

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