function [iseed,rand]=randf(iseed)
% Random number generation
% function [iseed,rand]=randf(iseed)
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
  mplier=16807; modlus=2147483647 ; mobymp=127773 ; momdmp=2836 ;
  nextn = iseed ;
  hvlue = fix(nextn / mobymp) ; lvlue = mod(nextn, mobymp) ;
  testv = mplier*lvlue - momdmp*hvlue ;
  if testv > 0
    nextn = testv ;
  else
    nextn = testv + modlus ;
  end
  rand = nextn/modlus ;
  iseed= nextn ;