www.pudn.com > NumericalComputingwithMatlabCode.zip > 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 
% See also RAND 
 
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); 
   mask = bitxor(bitshift(jhi,32,52),jlo); 
   [f,e] = log2(x); 
   x = pow2(bitxor(pow2(f,53),mask),e-53);