www.pudn.com > NumericalComputing.rar > randtx.m, change:2003-09-07,size:2600b

```function U = randtx(arg1,arg2)
% RANDTX  Text book version of RAND
% Uniformly distributed random numbers
% This M-file exactly reproduces the numerical
% behavior of the builtin RAND function.
% Usage:
%    randtx                  1-by-1
%    randtx(n)               n-by-n
%    randtx(m,n)             m-by-n
%    s = randtx('state')     Get state
%    randtx('state',j)       Set state determined by j
%    randtx('state',s)       Restore state

persistent z b i j ulp

if isempty(z)

% Initialize the state.

j = 2^31;
z = randsetup(32,j);
b = 0;
i = 0;
ulp = pow2(1,-53);
end

if nargin > 0 & isequal(arg1,'state')
if nargin == 1

% Get the state.

U = [z; b; pow2(i,-52); pow2(j,-52)];

elseif length(arg2) == 35

% Restore the state.

z = U(1:32);
b = U(33);
i = pow2(U(34),52);
j = pow2(U(35),52);

else

% Set the state determined by j.

j = arg2;
if j == 0, j = 2^31; end
z = randsetup(32,j);
b = 0;
i = 0;

end
return
end

% Determine size of output.

if nargin == 0
m = 1; n = 1;
elseif nargin == 1
m = arg1; n = arg1;
else
m = arg1; n = arg2;
end

% The MATLAB uniform random number generator.

U = zeros(m,n);
for k = 1:m*n
x = z(mod(i+20,32)+1) - z(mod(i+5,32)+1) - b;
if x < 0
x = x + 1;
b = ulp;
else
b = 0;
end
z(i+1) = x;
i = i+1;
if i == 32, i = 0; end
[x,j] = randbits(x,j);
U(k) = x;
end

% ----------------------------------------------------

function z = randsetup(n,j)
% RANDSETUP Generate n floats in [0,1] bit by bit.
% z = RANDSETUP(n,j)
% Seed is j.

z = zeros(n,1);
for k = 1:n
x = 0;
for d = 1:53
j = randint(j);
x = 2*x + bitand(bitshift(j,-19),1);
end
z(k) = pow2(x,-53);
end

% ----------------------------------------------------

function j = randint(j)
% RANDINT  Shift register random integer generator.
% The statement
%    j = randint(j)
% generates another random integer, 0 <= j <= 2^32-1

j = bitxor(j,bitshift(j,13,32));
j = bitxor(j,bitshift(j,-17,32));
j = bitxor(j,bitshift(j,5,32));

% ----------------------------------------------------

function [x,jhi] = randbits(x,jlo)
% RANDBITS  XOR random integer with floating fraction
%    [x,j] = randbits(x,j)

jhi = randint(jlo);