function r = randssp ( p, q )

%*****************************************************************************80
%
%% randssp() evaluates a multiplicative congruential uniform random number generator.
%
%  Discussion:
%
%    Based on the parameters used by IBM's Scientific Subroutine Package.
%    The statement
%      r = randssp
%    generates a single uniformly distributed random number.
%
%    The statement
%      r = randssp(m,n)
%    generates an m-by-n random matrix.
%
%    The statement
%      clear randssp
%    will cause the generator to reinitialize itself.
%
%    This function uses the "bad" generator parameters that IBM
%    used in several libraries in the 1960's.  There is a strong
%    serial correlation between three consecutive values.
%
%    See also RANDGUI, RANDMCG.
%
%  Licensing:
%
%    Copyright 2014 The MathWorks, Inc.
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625,
%    ebook: https://www.mathworks.com/moler/chapters.html
%
persistent m a c x
if isempty(x)
   m = 2^31;
   a = 2^16+3;
   c = 0;
   x = 123456789;
end

if nargin < 1, p = 1; end
if nargin < 2, q = p; end
r = zeros(p,q);
for k = 1:p*q
   x = rem(a*x + c, m);
   r(k) = x/m;
end
