Code covered by the BSD License

# Variable Precision Integer Arithmetic

### John D'Errico (view profile)

19 Jan 2009 (Updated )

Arithmetic with integers of fully arbitrary size. Arrays and vectors of vpi numbers are supported.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

factorialfactors(n)
```function [facs,freps] = factorialfactors(n)
% factorialfactors: efficient computation of the prime factors of factorial(n)
% usage: [facs,freps] = factorialfactors(n)
%
% A number sieve is used to generate the list of factors
%
% Arguments: (input)
%  n     - scalar, non-negative integer
%
% Arguments: (output)
%  facs  - row vector, a list of the prime factors of
%       factorial(n), sorted in increasing orsder.
%
%  freps - row vector, containing the number of occurrences
%       of each listed prime factor in the resulting
%       factorial.
%
%
% Example usage:
% % We could get the factors of factorial(100) using
% % factor(factorial(vpi(50))). But this is inefficient
% % and unnecessary.
%
% [facs,freps] = factorialfactors(50);
% [facs;freps]
% % ans =
% %    2  3  5  7 11 13 17 19 23 29 31 37 41 43 47
% %   47 22 12  8  4  3  2  2  2  1  1  1  1  1  1
%
% % Are these factors correct? The test is easy.
%
% factorial(vpi(50))
% % ans =
% %    30414093201713378043612608166064768844377641568960512000000000000
%
% prod(vpi(facs).^freps)
% % ans =
% %    30414093201713378043612608166064768844377641568960512000000000000
%
%
% % A second test for a bit larger argument
% tic
% [facs,freps] = factorialfactors(10000000);
% toc
%
% % Elapsed time is 2.277722 seconds.
%
% % Suppose we had tried to compute the actual
% % factorial itself, how many digits would it
% % (Answer: over 65 million decimal digits!)
%
% format long g
% sum(log10(ff).*fc)
%
% % ans =
% %           65657059.0800584
%
%

% check that n is a scalar, and a nonegative integer.
if (nargin ~= 1)
help factorialfactors
return
elseif isempty(n)
facs = [];
freps = [];
return
elseif (numel(n) ~= 1) || ~isnumeric(n) || (n<0) || (n ~= round(n))
error('FACTORIALFACTORS:improperargument', ...
'factorialfactors requires a scalar, nonnegative integer argument')
elseif (n==0) || (n==1)
% factorial(0) = factorial(1) = 1, so there
% are no prime factors.
facs = [];
freps = [];
return
end

% we have already eliminated n == 0 or 1. So all
% other cases have at least one prime factor.
facs = primes(n);
np = length(facs);
freps = zeros(1,np);

% start the sieving operation
for i = 1:np
p_i = facs(i);

% how many times will p_i^k appear as a factor?
ppower = p_i;
while ppower <= n
% remove factors of p_i from the product
ppreps = floor(n/ppower);
freps(i) = freps(i) + ppreps;

% incrementally multiply ppower
ppower = ppower*p_i;
end

end

```