Hump-day challenger - Recursion

In honor of John D'Errico (and some recent posts about recursion in MATLAB), I bring a recursion challenger.

Here is the challenge: Create a function which uses recursion to find the index location of one number in a vector of unique numbers. Your function should take a vector, and use recursion to find the index location of a value in the vector. The function should take both of the below values as arguments. The function should not call any built-in MATLAB set functions, including:

ISMEMBER, INTERSECT, SETDIFF, UNIQUE, UNION, SETXOR.  
Also, no calls to FIND or toolbox functions!

    M = randperm(10^4);  % The vector of unique numbers.
    V = ceil(rand*10^4) + 1;  % The value to find.

As you can see, the value V might not be in M. In this case your function should return an empty array. Remember, you must use recursion to do the work! And no fair increasing the recursion limit beyond 500!

Here is one solution to the problem, lets see how others do it!

Recursion Solution

[Jan]

Any further ideas? I thought recursion is an interesting field. I admit I avoid recursion in general, but obviously most people avoid it even better.

[Matt Fig]

I have another idea, but I don't have time to display it right now. Hope to get to it by Monday.

I agree that recursion in general is awful, and unnecessary. However, it is a good exercise to be able to think about recursive algorithms. So I also agree with John D'Errico's recommendation that noobs should definitely try to tackle a few recursive problems. Thanks Jan.

[Jan]

@Matt: I'm really curious about your idea, because I ran out of ideas completely. I'm surprised by the low traffic in this thread. I expected more comments or discussions...

[Jan]

Bump. I'm very curious about further recursive ideas.

[Walter Roberson]

I haven[t read any of the answers yet. The ideas that came to mind when I originally read the question should still work, but would probably have a limit of 2^500 elements in the list before overflowing the recursion. Are we excused from having to handle more than 2^53 elements?

[Matt Fig]

Walter, I wrote the rules as carefully as I could, with you in mind! As long as you can ddo 10^4 without hitting the limit, please post your solution!

[Matt Fig]

Thanks to Kenneth and Jan for their submissions!

Answer 1

[Kenneth Eaton]

My first attempt:

function index = find_index(M,V)

    switch numel(M)
      case 0
        index = [];
      case 1
        if M == V
          index = 1;
        else
          index = [];
        end
      otherwise
        nHalf = ceil(numel(M)/2);
        index = find_index(M(1:nHalf),V);
        if isempty(index)
          index = nHalf+find_index(M((nHalf+1):end),V);
        end
    end
end

[Matt Fig]

Nice, Kenneth. Different than the one I posted, which is good!

Answer 2

[Jan]

Kenneth's solution can be slightly modified to save memory:

function index = find_index(M, V, low, high)
if nargin == 2
  low = 1;
  high = length(M);
end

Len = high - low + 1;
switch Len
  case 0
    index = [];
  case 1
    if M(high) == V  % or M(low)
      index = high;
    else
      index = [];
    end
  otherwise
    nHalf = ceil(Len / 2);
    index = find_index(M, V, 1, nHalf);
    if isempty(index)
      index = find_index(M, V, nHalf + 1, high);
    end
  end
end

This cannot be called a new solution, but it was too ugly to post this without formatting in a comment. The underlying idea helps to avoid a common pitfall in recursive programming: wild data copy.

Answer 3

[Jan]

I admit one can discuss if a nested R(R(R(R(...)))) call is a "recursion" in the strict formal definition. It seems to satisfy the definition on Wikipedia. And it does not collide with Matlab RecursionLimit.

Unfortunately this program does not run in modern Matlab versions, which restrict the levels of nested parenthesis to 32, but in Matlab 6.5 this runs:

function Index = FindIndex   % Main function ---
N = 1e4;
M = randperm(N);
V = ceil(rand * N) + 1;

K = (V == M);
Data = [1, N+1, K];  % Need a single data vector

if any(K)
  s1(2:2:2*N) = '(';  % No REPMAT !
  s1(1:2:2*N) = 'R';
  s2(1:N)     = ')';
  Data = eval([s1, 'Data', s2]);
  Index = Data(2);
else
  Index = [];
end 

function Data = R(Data)  % Recursive subfunction: ---
if Data(1)  % Not found before
   Data(2) = Data(2) - 1;
   if Data(2 + Data(2))
      Data(1) = 0;
   end
end

Answer 4

[Jan]

If SPRINTF, SYSTEM and SAVE are accepted:

function R(M, V)
N = length(M);
if N == 0
  N = [];
  save('Result.mat', 'N');
else
  if M(end) == V
    save('Result.mat', 'N');
  else
    Cmd = ['R([', sprintf('%g ', M(1:end-1)), '], ', ...
          sprintf('%g', V), ')'];
    system(['matlab -r "', Cmd, '"']);
  end
end
quit;

Again a discussion is possible, if this is a valid recursion. But if you look inside the stack trace, even a standard call of an M-function starts a new instance of m_interpreter.dll -> _inInterPCode.

Answer 5

[Jan]

Because the vertical recursion would reach the RecursionLimit, we can split the numbers horizontally:

function Index = R2(M, V)
M(rem(M, 10) ~= rem(V, 10)) = -1;
match = (M >= 0);
if any(match)
  if sum(match) == 1 && V < 10
     K     = 1:length(M);
     Index = K(match);
  else
     M(match) = round(M(match) / 10);
     Index    = R2(M, round(V / 10));
  end
else  % No matching element:
  Index = [];
end

This function compares the rightmost digit and sets not matching values to -1. The depth of recursion is limited to LOG10(max(M)). An equivalent approach prints the values to a CHAR matrix at first and crops the righmost caracter in each recursion.

Answer 6

[Jan]

Flush the recursion at a certain level of recursion (< RecursionLimit, here: 100) and restart it:

function Index = FindIndex(M, V)
Index = R3(M, V, 0, 0);

% --- Recursive function: ---
function [Index, Found] = R3(M, V, Index, Level)

Index = Index + 1;
if Index > length(M)         % Last element reached:
   Index = [];
   Found = 1;   
elseif M(Index) == V         % Matching element found:
   Found = 1;
elseif Level < 100           % Call recursion
   [Index, Found] = R3(M, V, Index, Level + 1);
   if Found == 0             % Recursive call not successful:
      if Level == 0          % Restart recursion in base level only:
         [Index, Found] = R3(M, V, Index, 0);
      end
   end
else                         % Recursion limit reached:
   Found = 0;
end

Answer 7

[Matt Fig]

O.k., here is my other idea. It first sorts the array then uses a much more efficient recursive method than the algorithm I posted to start this thing off. Note that if it is objected that I used SORT, I could have written a recursive sort algorithm too (perhaps the subject of another challenger)! (Does this get me off the hook?) In my timings, given the inputs I included in the challenge description, this is about 70 times faster.

function idx = recurse_find_s(A,V)
% Returns the location of value V in a vector of unique values.
% A is the array, V the value to find.

[A,I] = sort(A);  
idx = I(find_array_sub(A,V,1,length(A)));

function idx = find_array_sub(A,V,L,R)
% Subfunction for recurse_find_s.  Recursion ahead.
if L>R
    idx = [];
else
    M = floor((L+R) /2);
    if V==A(M)
        idx = M;
    elseif V < A(M);  % Recursive call.
        idx = find_array_sub(A,V,L,M-1);
    else
        idx = find_array_sub(A,V,M+1,R);
    end
end

[Jan]

Very nice and "actually obvious" (this means: "obvioulsy not obvious"). Creating a binary search tree reduces the recursion depth substantially. +1

Answer 8

[Matt Fig]

O.k., here is something I do not understand. In trying to make the original algorithm I published (reproduced below) more efficient, I thought, "Why not use a nested function instead of a subfunction, that way we wouldn't have to pass around copies of the array." However, this is much slower! Can anyone explain why?

The original: .

.

function idx = recurse_find(A,V)
% Returns the location of value V in a vector of unique values.
% A is the array, V the value to find.

idx = find_array_sub(A,V,1,length(A));

function idx = find_array_sub(A,V,L,R)
% Subfunction for recurse_find.  Recursion ahead.
if L>R
    idx = [];
else
    M = floor((L+R) /2);
    if V==A(M)
        idx = M;
    else
        idxat = find_array_sub(A,V,L,M-1);  % Recursive call.
        if ~isempty(idxat)
            idx = idxat;
        else
            idx = find_array_sub(A,V,M+1,R) ;  % Recursive call.
        end

      end
  end

.

And the nested version of the same exact algorithm:

.

function idxn = recurse_find_n(A,V)
% Returns the location of value V in a vector of unique values.
% A is the array, V the value to find.

idxn = find_array_nest(1,length(A));

      function idx = find_array_nest(L,R)
      % Nested function, here is the recursion.
          if L>R
              idx = [];
          else
              M = floor((L+R) /2);
              if V==A(M)
                  idx = M;
              else
                  idxat = find_array_nest(L,M-1);
                  if ~isempty(idxat)
                      idx = idxat;
                  else
                      idx = find_array_nest(M+1,R) ;
                  end
              end
          end
      end
  end

[Jan]

Strange effect. I'd accuse the JIT to be too lazy inside a nested function. Passing copies of the arrays has the benefit, that there are no side effects.