memory overflow with double factorial function

I am trying to create a recursive formula to calculate the double factorial of a number, n. The code I have come up with currently is as follows:

function [DFact] = DFactorialRec(n) 
   DFact = n*DFactorialRec(n-2) 
end

This clearly does not work as it will recursively continue for negative values of n. Please could someone suggest an alteration to stop the memory overflow error?

================================================================

Out of memory. The likely cause is an infinite recursion within the program.

Error in DFactorialRec (line 2)

DFact = n*DFactorialRec(n-2)

================================================================

Thanks!

[David Goodmanson]

Hi Harry,

this is one of the the standard issues with recursion: when to stop. You just need to insert an if statement so that if n = 1 (having started with odd n originally) or n = 0 (having started with even n originally) then the output of DFactorialRec equals 1, and otherwise you do the recursive call that you already have.

[Harry Rogers]

hello David, thanks for the reply. I have updated my code as followed however I am now recieving a new error.

function [DFact] = DFactorialRec(n) 
if n==0
    disp("n is equal to 1")
else if n==1
        disp("n is equal to 1")
    else
   DFact = n*DFactorialRec(n-2) 
    end
end
end

the new error is as follows

Output argument "DFact" (and maybe others) not assigned during call to "DFactorialRec".

Error in DFactorialRec (line 7)

DFact = n*DFactorialRec(n-2)

any help would be greatly appreciated. Thanks

[Walter Roberson]

function [DFact] = DFactorialRec(n) 
if n==0
    disp("n is equal to 1")
    DFact = 1;
else if n==1
        disp("n is equal to 1")
        DFact = 1;
    else
   DFact = n*DFactorialRec(n-2);
    end
end
end

However I would recommend you rewrite using elseif instead of else if

[Harry Rogers]

thank you for the correction however this does not solve the error.. do you have any advice regarding how to make this work?

[Walter Roberson]

The code I posted worked for me provided that n is a scalar. If it is a non-scalar then you have the problem that your if is comparing a non-scalar to a scalar, and that comparison might be true for some elements but false for others...

Answer 1

[John D'Errico]

Recursive formulas like this always seem useful, but they are terrible in terms of coding, in terms of efficiency. You don't want to write it that way! Yes, the code you wrote looks pretty. That it will not work as you wrote it is a problem of course.

Think about it. For every recursive call, you need to allocate a new workspace, you incur function call overhead. This is a bad thing, and done for absolutely no good reason. Instead, a simple loop is all you need. Start at the bottom.

But better even yet is to not bother with an explicit loop. You can write this function in not much more than one line of code. Just use prod. The rest of what I did was mostly to make it friendly.

function DFact = DFactorial(n)
  % Double factorial function, n!!
  %  https://en.wikipedia.org/wiki/Double_factorial 
  %
  % n!! = 1 for both n == 1 and n == 0.
  if isempty(n) || (numel(n) > 1) || (n < 0) || (mod(n,1) ~= 0)
    error('The sky is falling. n must be scalar, non-negative, integer.')
  end
  
  DFact = 1;
  if n > 1
    start = 1 + mod(n + 1,2); % caters for either parity of n
    DFact = prod(start:2:n);
  end
end

Will this code fail for large values of n? Of course. It is easy to overwhelm the dynamic range of a double. But that will happen for any code that uses double precision. An exact integer will not be produced for n > 29 when done in double precision, and inf will result above 300. But again, that will happen in any case.

It is also easy to write a version that handles vector or array arguments for n.

[Harry Rogers]

Thats brilliant thank you. You've given me a deeper understanding about different operators that can be used within matlab. I'm faily new to using it so thanks for the advice.

[Nate Crummett]

Vectorized code from above:

function DFact = DFactorial(n)
% Double factorial function, n!!
% https://en.wikipedia.org/wiki/Double_factorial
% 
% n!! = 1 for both n == 1 and n == 0
% Inputs must be vectors of positive, real integers
%
% Designed by John D'Errico on 9 Oct 2020, taken from the MATLAB forum
% Vectorized by Robert Nate Crummett 30 July 2022

DFact = zeros(size(n));
idx_zeros = find(n == 0);
idx_ones = find(n == 1);
DFact([idx_ones, idx_zeros]) = 1;

n_end = n(~DFact);
start = 1 + mod(n_end + 1,2); % caters for either parity of n
prod_results = zeros(size(n_end));

for i = 1:length(n_end)
    prod_results(i) = prod(start(i):2:n_end(i));
end

DFact(~DFact) = prod_results;

end

Answer 2

[SteveC]

The double factorial of an integer -1 and larger is defined as in

https://mathworld.wolfram.com/DoubleFactorial.html

and generalized in the Pochhammer symbol

https://mathworld.wolfram.com/PochhammerSymbol.html

It appears in electromagnetic expansions and in special functions.

if the argument is an integer -1 and larger, it is written in Matlab as

fDF = @(bb) prod(bb:(-1):1)

For a vector example

>> arrayfun(fDF,-1:5)

1 1 1 2 6 24 120

If the application manages the validity of the argument, a single line without parity check is OK to define the function.

For big args, a Stirling form...

[SteveC]

Pardon the typo: (-2) not (-1) in the iterator.

>> fDF = @(bb) prod(bb:(-2):1)

fDF = @(bb)prod(bb:(-2):1)

> arrayfun(fDF,-1:10)

Columns 1 through 7

1 1 1 2 3 8 15

Columns 8 through 12

48 105 384 945 3840

>>