Seeding the pseudo random number generator

[Marc van Woerkom]

Am Wed,  9 Jun 2004 02:10:32 CEST hat

>Have a look to the Primes numbers generator from Joe's 
>example at 
>(http://www.erlang.org/examples/examples-2.0.html).
>
>The seed function is:
>
>new_seed() ->
>          {_,_,X} = erlang:now(),
>          {H,M,S} = time(),
>          H1 = H * X rem 32767,
>          M1 = M * X rem 32767,
>          S1 = S * X rem 32767,
>          put(random_seed, {H1,M1,S1}).
>
>This function generates different numbers each time when 
>it's called.

The ICFP contest featured such a modular arithmetic based 
pseudo random number generator as well.
See the part number theory in the task description.

   http://www.cis.upenn.edu/proj/plclub/contest/task.php

For me, struggling to get away from imperative paradigm to 
think in means functional programming, it posed some 
questions:

Functional programs are supposed to be side effect free, 
so how the heck do I ensure that I get a new number each 
time I called, that stuff is designed to give the same 
answer everytime.

Ok, Joe's examples does it by cheating, meaning calling a 
function with side effects (system time or such).

My only way out was to put in external information via an 
argument, which is used to carry the present seed series 
element from one call to another.

How to implement that series of Seeds S_i and pseudo 
random numbers X_i = F(S_{i+4})?

I needed to check against 100 sample randoms and thus came 
up with a first version A) that constructed a list 
[0,..,99] and then transformed that into a list, using 
list comprehension:

     [ begin {X, S} = randomint(N, I), io:format("~w, ", 
[X]) end || I <- lists:seq(0, 99)]

I looked that generator style up in the programming 
examples, interesting syntax, but it doesn't seem very 
space/run time efficient. 

The second attempt B) is based on the infinite lists 
example, making use of an lazy evaluation technique.
I think the list comrpehension syntax is more beautiful, 
but this beast is more efficient.

Any comments are welcome. I would love to know how you 
would code this.

Regards,
Marc


%
%  ICFP 2004
%  number_theory.erl
%

-module(number_theory).
-author('marc.vanwoerkom@REDACTED').

-include("definitions.hrl").

-export([s_from/1, 
	 x_from/1,
	 test_a/1, 
	 test_b/1]).

% random number generator (so that everyone uses the same 
one)

%
%  A) List comprehension version:
%

% this works faster, if one remembers the old seed S_old
old_randomint(N, S_old) ->
     S = (S_old * 22695477 + 1) rem (1 bsl 30),
     X = (S div 65536) rem 16384,
     {X rem N, S}.


% i-th random number from  {0, .., N-1}
randomint(N, 0) ->
     S0 = ?S0,
     {_, S1} = old_randomint(N, S0),
     {_, S2} = old_randomint(N, S1),
     {_, S3} = old_randomint(N, S2),
     old_randomint(N, S3);

randomint(N, I) ->
     {_, S_old} = randomint(N, I - 1),
     old_randomint(N, S_old).


% test 
test_a(N) ->
     [ begin {X, S} = randomint(N, I), io:format("~w, ", 
[X]) end || I <- lists:seq(0, 99)],
     ok.


%
%  B) lazy evaluation version:
%

% S series
s_from(N) ->
     s_from(N, ?S0).

% returns a fun that returns the pair {S_i, fun that 
returns next pair}
% takes numbers N and S_{i-1} as arguments
s_from(N, S) ->
     S1 = (S * 22695477 + 1) rem (1 bsl 30),  % bound S_i 
growth
     fun() ->
	    {S1, s_from(N, S1)}
     end.


% X series
x_from(N) ->
     FS1 = s_from(N),   % S generator fun FS1 returns S_1 
and next generator
     {_, FS2} = FS1(),
     {_, FS3} = FS2(),
     {_, FS4} = FS3(),
     x_from(N, FS4).    % X_0 will be calculated from S_4 
etc

% returns a fun that returns the pair {X_i, fun that 
returns next pair}
% takes N and old S generator fun FS as arguments
x_from(N, FS) ->
     {S, FS1} = FS(),
     X = (S div 65536) rem 16384,
     fun() ->
	    {X rem N, x_from(N, FS1)}
     end.


% test 100 nums
test_b(N) -> test_b(N, x_from(N), 0, 100).

test_b(N, _, _, 0) ->
     io:format(" ok~n");

test_b(N, F, I, J) ->
     {X, F1} = F(),
     io:format("X~w = ~w, ", [I,X]),
     test_b(N, F1, I+1, J-1).