[runifresult,t,n]=runif(t,n);

function [runifresult,t,n]=runif(t,n);
runifresult=[];
persistent dummy firstCall floatn i j nold runif ; if isempty(firstCall),firstCall=1;end; 

if isempty(dummy), dummy=0; end;
if isempty(floatn), floatn=0; end;
if isempty(runifresult), runifresult=0; end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(nold), nold=0; end;
%***BEGIN PROLOGUE  RUNIF
%***PURPOSE  Generate a uniformly distributed random number.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  L6A21
%***TYPE      SINGLE PRECISION (RUNIF-S)
%***KEYWORDS  FNLIB, RANDOM NUMBER, SPECIAL FUNCTIONS, UNIFORM
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% This random number generator is portable among a wide variety of
% computers.  It generates a random number between 0.0 and 1.0 accord-
% ing to the algorithm presented by Bays and Durham (TOMS, 2, 59,
% 1976).  The motivation for using this scheme, which resembles the
% Maclaren-Marsaglia method, is to greatly increase the period of the
% random sequence.  If the period of the basic generator (RAND) is P,
% then the expected mean period of the sequence generated by RUNIF is
% given by   new mean P = SQRT (PI*FACTORIAL(N)/(8*P)),
% where FACTORIAL(N) must be much greater than P in this asymptotic
% formula.  Generally, N should be around 32 if P=4.0E6 as for RAND.
%
%             Input Argument --
% N      ABS(N) is the number of random numbers in an auxiliary table.
%        Note though that ABS(N)+1 is the number of items in array T.
%        If N is positive and differs from its value in the previous
%        invocation, then the table is initialized for the new value of
%        N.  If N is negative, ABS(N) is the number of items in an
%        auxiliary table, but the tables are now assumed already to
%        be initialized.  This option enables the user to savemlv the
%        table T at the end of a long computer run and to restart with
%        the same sequence.  Normally, RUNIF would be called at most
%        once with negative N.  Subsequent invocations would have N
%        positive and of the correct magnitude.
%
%             Input and Output Argument  --
% T      an array of ABS(N)+1 random numbers from a previous invocation
%        of RUNIF.  Whenever N is positive and differs from the old
%        N, the table is initialized.  The first ABS(N) numbers are the
%        table discussed in the reference, and the N+1 -st value is Y.
%        This array may be saved in order to restart a sequence.
%
%             Output Value --
% RUNIF  a random number between 0.0 and 1.0.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  RAND
%***REVISION HISTORY  (YYMMDD)
%   770401  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   910819  Added EXTERNAL statement for RAND due to problem on IBM
%           RS 6000.  (WRB)
%***end PROLOGUE  RUNIF
t_shape=size(t);t=reshape(t,1,[]);
if firstCall,   nold=[-1];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  RUNIF
if( n~=nold )
%
nold = fix(abs(n));
floatn = nold;
if( n<0 )
[ dummy ,t(nold+1)]=rand(t(nold+1));
end;
if( n>=0 )
%
for i = 1 : nold;
[t(i) ]=rand(0.);
end; i = fix(nold+1);
[t(nold+1) ]=rand(0.);
end;
end;
%
j = fix(t(nold+1).*floatn + 1.);
t(nold+1) = t(j);
runifresult = t(j);
[t(j) ]=rand(0.);
%
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',t); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',n); evalin('caller',[inputname(2),'=FUntemp;']); end
end
%DECK RWUPDT