%%
%% Simple random number generator used for tests in order to generate the
%% same sequence of random numbers under Octave and Matlab
%%

function y = my_rand(varargin)

global seed xn;
a = int32(16807);
c = int32(1013904223);
m = int32(2^31-1);

disp(sprintf('State at begin: a=%d, c=%d, m=%d, seed=%d, xn=%d', a, c, m, seed, xn));

% Determine the calling mode
if nargin == 0 % Generate 1 random number
    mode = 'generate';
    s = [1 1];
else
    if ischar(varargin{1}) % Set or get the seed
        if ~strcmp(varargin{1}, 'state')
            error('I''m too stupid, use ''state''');
        end
        if nargin == 1 % return the seed
            mode = 'getseed';
        else % Reset seed
            mode = 'setseed';
            if numel(varargin{2}) ~= 1
                error('One seed at a time');
            end
        end
    else
        mode = 'generate';
        if nargin == 1 % all size in one arg
            if isscalar(varargin{1})
                s = [varargin{1}, varargin{1}];
            else
                s = varargin{1};
            end
        else % size in different args
            s = cell2mat(varargin);
        end
    end;
end;

switch mode
    case 'getseed' % Return the seed
        y = seed;
    case 'setseed' % Reset seed
        seed = int32(varargin{2});
        %n = int32(0);
        xn = int32(rem(a*seed, m));
    case 'generate'
        disp(sprintf('Gonna generate %d elements', prod(s)));
        y = int32(zeros(s));
        y(1) = xn;
        for i=2:numel(y)
            y(i) = int32(rem(a*y(i-1), m));
        end
        xn = int32(rem(a*y(end), m));
        %y = reshape(y, s);
    otherwise
        error('Too dumb');
end;

disp(sprintf('State at end: a=%d, c=%d, m=%d, seed=%d, xn=%d', a, c, m, seed, xn));