myfactor(number)

function [output] = myfactor(number)

%MYFACTOR Prime factors of large integers 
%
% PURPOSE:
%
% This function determines all the factors of the input number.
% The program also indicates if the number is prime.
%
% For numbers up to 10^14, an accurate result is guaranteed.
%
% For numbers greater than 10^14, it may work. For certain cases (which will 
% be signalled by the program) the result cannot be guaranteed without 
% changing the code to use more CPU time and computer memory.
% If the program is unhappy about the result for certain numbers > 10^14,
% it will say so.
%
% Procedure for numbers > 10^14: First, run the program in its present form.
% In most cases, the answer will be straightforward. 
% If the program announces a problem, you must consider increasing the
% variable (limit_parameter) in line 52 of the program from its current value
% of 10^7, up to the square root of the largest factor in the answer, but this 
% can make significant demands on computer resources.
%
% SYNOPSIS: myfactor(number)
%
% EXAMPLE:  myfactor(2^32+1) 
%           ans = [1 641 6700417]
%

if numel(number)~=1, error('MATLAB:factor:NonScalarInput','input must be a scalar.'); end
if (number < 1) || (floor(number) ~= number)
  error('MATLAB:factor:InputNotPosInt', 'input must be a positive integer.'); 
end

if number <= 2^32
    disp(' '); disp('This input is in the domain [0 2^32] and is solved directly using MATLAB''s built-in function factor.m')
    output = factor(number);
    if length(output == 1) 
       if output > 1, output = [1 output]; end
    end
    return
end

if number > 10^14,
   disp(' '),
   disp('For inputs below 10^14, this program guarantees a result.'),
   disp('For inputs above 10^14, it makes an effort.') 
   disp(' '),
end

limit_parameter = 10^7;
divisor = [1 primes(min(floor(sqrt(number)),limit_parameter))];
range =  2 * length(divisor);
test_factor(1) = number; test_factor(2) = 1; 
quotient(1) = number; 
factor_count = 2; prime_count = 2; 

for i = 2 : range
    if   rem(quotient(factor_count-1), divisor(prime_count)) == 0, 
         test_factor(factor_count+1) = divisor(prime_count);
         quotient(factor_count) = quotient(factor_count-1)/divisor(prime_count);
         factor_count = factor_count + 1;
    else prime_count = prime_count + 1;
         if prime_count > length(divisor), break; 
         end
    end
end

if     factor_count == 2
       output = [1 number];
elseif quotient(factor_count-1) > 1
       test_factor(factor_count+1) = quotient(factor_count-1);
       output = test_factor(3:factor_count+1);
else   output = test_factor(3:factor_count);         
end

if   min(output) == 1,
     if max(output) <= limit_parameter^2
      disp(' '),disp([int2str(number),' is a prime number.'])
     end
else output = [1 output]; 
end

if max(output) > limit_parameter^2
disp('In this case, the program is unable to test the largest factor shown in the answer.'),
disp('In order to do this, it would be necessary to increase the variable (limit_parameter)'),
disp('in line 52 of the program, presently set at 10^7, up to the square root of the largest'),
disp('factor, but this could have consequences in terms of CPU time and computer memory.'),
end

if max(output) > 10^9
output = uint64(output);
end