randword(alpha, wordSize, numWords)

function word = randword(alpha, wordSize, numWords)
%RANDWORD Random word from alphabet
%    WORD = RANDWORD(alpha, wordSize) returns a random word in the alphabet
%    'alpha', with the number of symbols 'wordSize'. The result WORD has
%    size 1 x wordSize.
%
%    WORD = RANDWORD(alpha, wordSize, numWords) returns 'numWords' random
%    words in the alphabet 'alpha', with the number of symbols 'wordSize'.
%    The result WORD has size numWords x wordSize.
%
%    Example:  
%    Generate 4 words of length 5 from the alphabet 'A', 'B', ..., 'Z'
%       word = randword('A':'Z', 5, 4)
%       word =
%       VPQQT
%       PFQPQ
%       OHPYM
%       XWAPM
%
%    The method of randword is to select a random word from the
%    lexicographically sorted list of all words on alpha of length
%    wordSize. With n = length(alpha) and w = wordSize, the steps are
%       1) Get a random number r on the interval 1:n^w.
%       2) Compute base n expansion of r.
%       3) Make word with symbols of 'alpha' indexed by base n digits of r.
%
%   Author: Tom Richardson
%           trichardson AT xrite DOT com
%   Date:   14 July 2009
%
%   Changed the routine to use multiple passes to ensure that values of
%   lexIndex stay below 2^31 - 1, thus avoiding roundoff error in
%   alphaSize.
%   Modified: Tom Richardson
%   Date:   21 January 2019
%       

    % Sort alpha and compute the size of alpha
    alpha = sort(alpha(:));
    alphaSize = length(alpha);
    
    % Compute one word if number of words unspecified
    if nargin < 3
        numWords = 1;
    end
   
    subwordSize = floor(31/log2(alphaSize));
    
    word = [];
    
    for z = 1:ceil(wordSize/subwordSize)
    
        % Compute random index into all words
        lexIndex = floor((alphaSize^subwordSize)*rand(numWords,1));

        % Find the index_th word in lexicographic order by selecting
        % the symbols of alpha indexed by the base alphaSize digits of lexIndex.
        % Compute base alphaSize digits of lexIndex, add 1 to index alpha.
        alphaIndex = 1 + mod(floor(bsxfun(@rdivide, lexIndex, alphaSize.^(0:subwordSize-1))), alphaSize);

        % Make word by indexing into alpha
        word = [word alpha(alphaIndex)];
    end
    
    word(:,wordSize+1:end) = [];
    
end