stirling.m

% compute the n-th order stirling numbers
% of second species for n>= 2
% 
% octave has a hard limit on recursion; 256 calls are the most
% you can make without raising an error
% set max_recursion_depth to your value
function x = stirling_recursive (n)
  if (n == 2)
    x =[1,1];
    return;
  endif
  % compute the previous order solution
  y = stirling_recursive(n-1);
  a = 2 : n-1;
  % S_{i-1} + i S_{i}
  temp =  y(1:end-1) + (y(2:end) .* a);
  x = [1, temp, 1];
endfunction

#            Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------------
   1  stirling_recursive    R        0.027      89.83           99
   5           binary .*             0.001       2.97           98
   3            binary -             0.001       2.28          294
   4                 end             0.001       1.93          196
   6            binary +             0.000       1.19           98
   2           binary ==             0.000       1.05           99
   7             profile             0.000       0.63            1
  10               false             0.000       0.05            1
   8              nargin             0.000       0.04            1
   9           binary !=             0.000       0.02            1
  11 __profiler_enable__             0.000       0.00            1

  
% compute the n-th order stirling numbers
%
function stir = stirling_iterative (n)
  % x :: [i32; n+1]
  x = zeros(1,n+1);

  % initialise stirling numbers of order 2
  x(2:3) = 1;

  if n == 2
    stir = x(2:3);
    return
  endif

  % x :: [i32; n+1]
  a = 0: n;

  for i = 2 : n-1
    % S_{i-1} + i S_{i}
    %soluzione cacca
    %x = shift(x,1) + x .* a;
    x = [0, x(1:end-1)] + x .* a;
  endfor

  stir = x(2:end);
endfunction

   #            Function Attr     Time (s)   Time (%)        Calls
------------------------------------------------------------------
   1  stirling_iterative             0.013      82.07            1
   7           binary .*             0.001       7.02           98
   2            binary +             0.001       3.86           99
   6                 end             0.000       3.07           99
   5            binary -             0.000       2.30           99
   8             profile             0.000       1.27            1
   3               zeros             0.000       0.15            1
  11               false             0.000       0.10            1
   9              nargin             0.000       0.07            1
  10           binary !=             0.000       0.04            1
   4           binary ==             0.000       0.03            1
  12 __profiler_enable__             0.000       0.00            1

soluzione cacca
   #           Function Attr     Time (s)   Time (%)        Calls
-----------------------------------------------------------------
  18             repmat             0.045      39.04           98
   6              shift             0.040      35.04           98
   1 stirling_iterative             0.014      12.32            1
  20                all             0.002       1.34          294
  15               size             0.001       1.21          294
  16           binary >             0.001       1.11          392
  17               find             0.001       1.06           98
   7             nargin             0.001       1.06          393
  19            isempty             0.001       0.84          392
   4          binary ==             0.001       0.82          491
  25            reshape             0.001       0.62           98
  10           binary <             0.001       0.58          294
  24               prod             0.001       0.56           98
  27          binary .*             0.001       0.55           98
  22                max             0.001       0.53           98
  23                any             0.001       0.45           98
   2           binary +             0.000       0.41          197
   9              numel             0.000       0.37          196
  13           prefix !             0.000       0.33          196
  26                rem             0.000       0.29           98