subfactorial(n)

function [r] = subfactorial(n)
%% subfactorial : computes the subfactorial of a number
%
% Usage :
%   r = subfactorial(n) calculates the subfactorial of a non negative integer n
%       using the inclusion-exculsion principle.
%   The subfactorial can be written as,
%       !n = n! sum( (-1)^k / k! ), for k = 0 to n
%       !n : subfactorial(n)
%       n! : factorial(n) = 1*2*3*...*n
%
% Remark :
%   This function is limited by largest floating point number that can be
%   defined, as given by the realmax. Since,
%       subfactorial(170) < realmax < subfactorial(171)
%   Maximum value of input (n) cannot exceed 170.
%   Higher value of (n) will return NaN.
%
% Example:
%   a = subfactorial(12);
%   returns -> a = 176214841
%
% See Also:
%   factorial

% Author: Yash Kochar ( [e-mail] )
% Last modified: 25th June 2009
%-------------------------------------------------------------------------------

% Check input :
if ~isscalar(n)
    error('subfactorial:NonScalarInput','Input must be a scalar.');
end
if (n < 0) || (floor(n) ~= n)
  error('subfactorial:NonNegativeIntegerOnly', ...
      'Input must be a non-negative integer.'); 
end
if n > 170
    r = NaN;
    return
end

% Calculate subfactorial :
a = fliplr( cumprod( [1, n: -1: 1]));
b = -2 *mod(0: n, 2) +1;
r = sum(a.*b);

% Alternate implementations :
% Also using inculsion-exculsion principle (slower)
% r = sum((fliplr(cumprod([1, n:-1:1]))).*((-1).^(0:n)))
%
% Direct computation (fastest, but not sure it if it is always correct)
% r = floor(factorial(n)/exp(1) +0.5);
%-------------------------------------------------------------------------------