% custom randomizer 
% must define seed
% x_n+1 = 8121*x_n + 28411, then (mod 134456)
% force dist between 0 and 1
% A. Sun 2/2005

pts=1000;
a=8121;
b=28411;
c=134456;


n=1:pts;
rand_vec=[seed];

for i=1:pts,
  r=mod((a*rand_vec(i)+b),c);
  rand_vec=[rand_vec r];
end

rand_vec = rand_vec / c;
% plot(rand_vec,'b.');