function [L,M,X,S] = logfactorial(N,varargin)
% LOGFACTORIAL computes log10 of the factorial.
%
% USAGE:
%   logfactorial(N)
%   logfactorial(N,method)
%   logfactorial(N,method,fmt)
%   logfactorial(N,'sum',fmt,lengthOfSequence)
%   [M,X,L,S] = logfactorial(N,...);
%
% INPUT:
%   N is a positive integer or an array of positive integers.
%
%   method is either 'gamma' (the default) or 'sum'. It selects the method
%   that is used for computation.
%
%   fmt is the format in which the mantissa in S is printed. This argument
%   is optional; default format is '%g'.
%
%   lengthOfSequence is the length of a vector that MATLAB is able to
%   assign. One can set this to Inf, but this entails the risk of running
%   out of memory. The default is 1E7. If MATLAB runs out of memory, choose
%   a smaller number. This option is only relevant when using the 'sum'
%   method.
%
% OUTPUT:
%   L is an array of doubles containing log10 of the factorials.
%   M is an array of doubles containing the mantissas of the results.
%   X is an array of integers containing the exponents of the results.
%   S is an array of cellstrings that represents the results as strings in
%   the form 'Me+X'.
%
%   The factorial of N is M .* 10.^X, or equivalently 10.^L. Consequently,
%   10.^logfactorial(N) produces the same result as MATLAB's standard
%   function factorial(N).
%
% EXAMPLE 1
%   logfactorial(1E6)
%   produces 5.5657e+006, which is log10(1E6!).
%
%   [L,M,X,S] = logfactorial(1E6)
%   produces
%   L = 5.5657e+006         % log10(N!)
%   M = 8.2639              % mantissa
%   X = 5565708             % exponent
%   S = '8.26393e+5565708'  % N! as a string
%
% EXAMPLE 2
%   [L,M,X,S] = logfactorial(21:25)
%   returns
%   L =  19.7083   21.0508   22.4125   23.7927   25.1906
%   M =   5.1091    1.1240    2.5852    6.2045    1.5511
%   X =  19        21        22        23        25
%   S =  '5.10909e+19' '1.124e+21' '2.5852e+22' '6.20448e+23' '1.55112e+25'
%
% REMARK: The 'gamma' method computes the log (wrt base 10) of the
% factorial by means of the gamma function. The 'sum' method computes the
% log of the factorial as sum(log10(i)), where i = 1 ... n. The
% exponential of the factorial is just the integer part of log10 of the
% factorial; the mantissa is 10 to the fractional part.
% The 'gamma' method is generally faster and more accurate, and should
% therefore always be used. The 'sum' method is provided here only for
% comparison, and because the algorithm is quite transparent.
%
% Author: Yvan Lengwiler (yvan.lengwiler@unibas.ch)
% Date: May 9, 2007
%
% Acknowledgement: Urs Schwarz and, in particular, John D'Errico (reviewers
% at MathCentral) for helpful comments.
%
% See also: FACTORIAL, LARGEFACTORIAL (MathCentral file id 14816)

% Compare the two methods with factorial.m for N = 1 ... 170 as follows:
%   dgamma = @(x) log10(factorial(x)) - logfactorial(x,'gamma');
%   dsum = @(x) log10(factorial(x)) - logfactorial(x,'sum');
%   hold all
%   plot(dgamma(1:170))
%   plot(dsum(1:170))
%   legend('gamma method','sum method')
% One can see the differences in precision already at these small values of
% N. The loss of precision gets worse for larger arguments.

    % --- check validity of arguments
    
    Nsort = N(:);   % This is not sorted yet, but will be sorted later. The
                    % same variable name is used to save memory.

    if any(fix(Nsort) ~= Nsort) || any(Nsort < 0) || ...
            ~isa(Nsort,'double') || ~isreal(Nsort)
        error('N must be a matrix of non-negative integers.')
    end
    
    if length(varargin) > 3
        error('Too many input arguments.')
    end
    
    % --- assign default arguments if necessary
    % (this trick is taken from Nabeel Azar, see
    % http://www.ee.columbia.edu/~marios/matlab/Programming%20Patterns%20Think%20Globally,%20Act%20Locally.pdf)
    
    defaultValues                 = {'gamma','%g',1E7};
    nonemptyIdx                   = ~cellfun('isempty',varargin);
    defaultValues(nonemptyIdx)    = varargin(nonemptyIdx);
    [method fmt lengthOfSequence] = deal(defaultValues{:});
    
    % --- perform computation
    
    switch lower(method)
        
        case 'gamma'	% ... using gammaln

            L = gammaln(N+1)/log(10);   % John D'Errico one line solution
                                        % of the problem
            
        case 'sum'      % ... using direct summation

            % We compute the smallest factorial first, ...
            [Nsort,map] = sort(Nsort);
            Nsort = [0;Nsort];
            % ... and then compute only the increments.
            L = arrayfun(@sumlogIncrement, 1:numel(N));
            L = cumsum(L);  % aggregate increments
            L(map) = L;     % put results back into original ordering
                            % ("un-sort")
            L = reshape(L,size(N));     % bring output into correct shape

        otherwise

            error(['Method ''', method, ''' is unknown.'])

    end

    if nargout > 1
        % 10^L is the factorial we are after, so ...
        X = fix(L);      % ... X is the exponential of the factorial ...
        M = 10.^(L-X);   % ... and M is the mantissa.
    end
    
    if nargout > 3
        S = arrayfun(@makeString, 1:numel(N));
        S = reshape(S,size(N));     % bring output into correct shape
    end

    % --- nested function -------------------------------------------------
    function s = sumlogIncrement(j)
    % sumlogIncrement(j) computes sum(log10(i)) for i = Nsort(j)+1 ...
    % Nsort(j+1); this is used for the 'sum' method only.
        lengthInterval = Nsort(j+1)-Nsort(j);
        a = Nsort(j);
        if lengthInterval <= lengthOfSequence
            s = sum(log10(a+(1:lengthInterval)));
        else
            % MATLAB sometimes runs out of memory when allocating long
            % vectors. It is assumed that MATLAB can always allocate a
            % vector of size 'lengthOfSequence' (chosen by user, or set to
            % 1E7 by default).
            seq = 1:lengthOfSequence;   % "short" sequence
            rounds = floor(lengthInterval/lengthOfSequence);  % # of rounds
            partialSum = zeros(1,rounds);
            for r = 1:rounds
                partialSum(r) = sum(log10(seq+(r-1)*lengthOfSequence+a));
            end
            s = sum(partialSum) + ...
                sum(log10(((rounds*lengthOfSequence+1):lengthInterval)+a));
        end
    end

    function str = makeString(j)
    % makeString(j) creates a string in the form 'Me+X'.
        str = cellstr(strcat(num2str(M(j),fmt), 'e+', int2str(X(j))));
    end

end
