以下是六个不同实现的比较,用于查找整数的因素：

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