root/library/utia_legacy/ticket_12/randun.m @ 819

function r=randun(m,n)
%RANDUN generate sample from Uniform distribution
%
% r = randun
% r = randun(m)
% r = randun(m,n)
%
% r       : generated sample
% SEED    : global variable SEED
% m       : rows
% n       : columns
%
% Design  : J. Andrysek
% Project : BadDyr
% Comments: used to generate exactly the same numbers in .m and .dll
% versions
%
% Patched by Ludvik Tesar to output matrix instead of scalar
 

global SEED;
if(isempty(SEED)) SEED = 1111111; end;
%SEED = randn('seed');

global RANDUN_COUNTER;
if(isempty(RANDUN_COUNTER)) RANDUN_COUNTER = 0; end;


if ~exist('m','var'); m = 1; end;
if ~exist('n','var'); n = 1; end;

r = zeros(m,n);

for i=1:(m*n);

% Ahrens linear congruential pseudo-random generator with A=663608941 B=0 M=2^32
%    A      = 663608941;
%    M      = 4294967296;
%    SEED   = soucinek(A,SEED);
%    r(i)   = SEED/M;
% replaced by:

% Park-Miller linear congruential pseudo-random generator with A=16807 B=0 M=2^31-1
%     (it spans all 2^31-1 numbers and has good statistical properties)
    A    = 16807;
    M    = 2147483647;
    SEED = mod(A*SEED,M);
    r(i) = SEED/M;
end;

%randn('seed',SEED);

RANDUN_COUNTER = RANDUN_COUNTER+m*n;

function r=soucinek(a,b);
z=2^16;
ha(1)=floor(a/z);
ha(2)=mod(a,z);

hb(1)=floor(b/z);
hb(2)=mod(b,z);

c=ha(2)*hb(2);
res(2)=mod(c,z);
res(1)=mod(floor(c/z)+hb(2)*ha(1)+hb(1)*ha(2),z);

r=res(1)*z+res(2);