Linear congruential generator

Linear congruential generator is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The linear congruential generator is a very simple example of a random number generator. All linear congruential generators use this formula:

• $r_{n+1}=a\times r_{n}+c{\pmod {m}}$

Where:

• $r_{0}$ is a seed.
• $r_{1}$, $r_{2}$, $r_{3}$, ..., are the random numbers.
• $a$, $c$, $m$ are constants.

If one chooses the values of $a$, $c$ and $m$ with care, then the generator produces a uniform distribution of integers from $0$ to $m-1$.

LCG numbers have poor quality. Anyone who knows $r_{n}$ can predict $r_{n+1}$, therefore LCG is not cryptographically secure. The LCG is still good enough for simple tasks like the Knuth shuffle. Among the benefits of the LCG, one can easily reproduce a sequence of numbers, from the same $r_{0}$. One can also reproduce such sequence with a different programming language, because the formula is so simple.

The task is to replicate two historic random number generators. One is the rand() function from BSD libc, and the other is the rand() function from the Microsoft C Runtime (MSCVRT.DLL). Each replica must yield the same sequence of integers as the original generator, when starting from the same seed.

(These formulas should have better wiki formatting.)

BSD formula

• ${\begin{array}{lcl}seed&=&1103515245\times seed+12345{\pmod {2^{31}}}\\rand&=&seed\end{array}}$

Microsoft formula

• ${\begin{array}{lcl}seed&=&214013\times seed+2531011{\pmod {2^{31}}}\\rand&=&seed\div 2^{16}\end{array}}$

(Here should be links to the sources of the formulas. Here should explain how the low bits have a short cycle.)

(Here should be some sample numbers for a few seeds, so that one who solves this task can check if one's numbers matches these samples.)

Ada

We first specify a generic package LCG:

<lang Ada>generic

  type Base_Type is mod <>;
  Multiplyer, Adder: Base_Type;
  Output_Divisor: Base_Type := 1;

package LCG is

  procedure Initialize(Seed: Base_Type);
  function Random return Base_Type;
  -- changes the state and outputs the result

end LCG;</lang>

Then we provide a generic implementation:

<lang Ada>package body LCG is

  State: Base_Type := Base_Type'First;

  procedure Initialize(Seed: Base_Type) is
  begin
     State := Seed;
  end Initialize;

  function Random return Base_Type is
  begin
     State := State * Multiplyer + Adder;
     return State / Output_Divisor;
  end Random;

end LCG;</lang>

Next, we define the MS- and BSD-instantiations of the generic package:

<lang Ada>with Ada.Text_IO, LCG;

procedure Run_LCGs is

  type M31 is mod 2**31;

  package BSD_Rand is new LCG(Base_Type => M31, Multiplyer => 1103515245,
                              Adder => 12345);

  package MS_Rand  is new LCG(Base_Type => M31, Multiplyer => 214013,
                              Adder => 2531011, Output_Divisor => 2**16);

begin

  for I in 1 .. 10 loop
     Ada.Text_IO.Put_Line(M31'Image(BSD_Rand.Random));
  end loop;
  for I in 1 .. 10 loop
      Ada.Text_IO.Put_Line(M31'Image(MS_Rand.Random));
  end loop;

end Run_LCGs;</lang>

Finally, we run the program, which generates the following output (note that the first ten lines are from the BSD generator, the next ten from the MS generator):

 12345
 1406932606
 654583775
 1449466924
 229283573
 1109335178
 1051550459
 1293799192
 794471793
 551188310
 38
 7719
 21238
 2437
 8855
 11797
 8365
 32285
 10450
 30612

C

In a pretended lib style, this code produces a rand() function depends on compiler macro: gcc -DMS_RAND uses MS style, otherwise it's BSD rand by default. <lang C>#include <stdio.h>

/* always assuming int is at least 32 bits */ int rand(); int rseed = 0;

inline void srand(int x) { rseed = x; }

1. ifndef MS_RAND
2. define RAND_MAX ((1U << 31) - 1)

inline int rand() { return rseed = (rseed * 1103515245 + 12345) & RAND_MAX; }

1. else /* MS rand */

1. define RAND_MAX_32 ((1U << 31) - 1)
2. define RAND_MAX ((1U << 15) - 1)

inline int rand() { return (rseed = (rseed * 214013 + 2531011) & RAND_MAX_32) >> 16; }

1. endif/* MS_RAND */

int main() { int i; printf("rand max is %d\n", RAND_MAX);

for (i = 0; i < 100; i++) printf("%d\n", rand());

return 0; }</lang>

Icon and Unicon

The following were written before the task was complete and may need adjustment. Both behave in basically the same way. The LCRNG maintains the state (seed) from round to round. The seed can be changed by specifing an argument. Otherwise it defaults. Also these must be tested for 32 bit overflow (see Discussion / Talk pages). <lang Icon> procedure rand_BSD(x) static seed,m initial {

  seed := 0    # default if no x given??
  m := 2^31-1
  }
  
  seed := \x
  seed := (1103515245 * seed + 12345) % m
  return seed

end

procedure rand_MS(x) static seed,m,d initial {

  seed := 0    # default if no x given??
  m := 2^31-1
  d := 2^16
  }

  seed := \x
  seed := (214013 * seed + 2531011) % m
  return seed/d

end</lang>

J

Solution: <lang j>lcg=: adverb define

0 m lcg y

'a c mod'=. m
}. (mod | c + a * ])^:(<y+1) x 

)

rand_bsd=: (1103515245 12345 , <.2^31) lcg rand_ms=: (2^16) <.@:%~ (214013 2531011 , <.2^31) lcg</lang> Example Use: <lang j> rand_bsd 10 12345 1406932606 654583775 1449466924 229283573 1109335178 1051550459 1293799192 794471793 551188310

  654583775 rand_bsd 4

1449466924 229283573 1109335178 1051550459

  rand_ms 10

38 7719 21238 2437 8855 11797 8365 32285 10450 30612</lang>

PARI/GP

Note that up to PARI/GP version 2.3.0, random() used a linear congruential generator. <lang parigp>BSDseed=Mod(1,1<<31); MSFTseed=Mod(1,1<<31); BSD()=BSDseed=1103515245*BSDseed+12345;lift(BSDseed); MSFT()=MSFTseed=214013*MSFTseed+2531011;lift(MSFTseed)%(1<<31);</lang>

PicoLisp

<lang PicoLisp>(zero *BsdSeed *MsSeed)

(de bsdRand ()

  (setq *BsdSeed
     (& (+ 12345 (* 1103515245 *BsdSeed)) `(dec (** 2 31))) ) )

(de msRand ()

  (>> 16
     (setq *MsSeed
        (& (+ 2531011 (* 214013 *MsSeed)) `(dec (** 2 31))) ) ) )</lang>

Output:

: (do 7 (printsp (bsdRand)))
12345 1406932606 654583775 1449466924 229283573 1109335178 1051550459 -> 1051550459

: (do 12 (printsp (msRand)))
38 7719 21238 2437 8855 11797 8365 32285 10450 30612 5853 28100 -> 28100

Ruby

You can create multiple instances of LCG::Berkeley or LCG::Microsoft. Each instance privately keeps the original seed in @seed, and the current state in @r. Each class resembles the core Random class, but with fewer features. The .new method takes a seed. The #rand method returns the next random number. The #seed method returns the original seed.

<lang ruby>module LCG

 module Common
   # The original seed of this generator.
   attr_reader :seed

   # Creates a linear congruential generator with the given _seed_.
   def initialize(seed)
     @seed = @r = seed
   end
 end

 # LCG::Berkeley generates 31-bit integers using the same formula
 # as BSD rand().
 class Berkeley
   include Common
   def rand
     @r = (1103515245 * @r + 12345) & 0x7fff_ffff
   end
 end

 # LCG::Microsoft generates 15-bit integers using the same formula
 # as rand() from the Microsoft C Runtime.
 class Microsoft
   include Common
   def rand
     @r = (214013 * @r + 2531011) & 0x7fff_ffff
     @r >> 16
   end
 end

end</lang>

The next example sets the seed to 1, and prints the first 5 random numbers.

<lang ruby>lcg = LCG::Berkeley.new(1) p (1..5).map {lcg.rand}

1. prints [1103527590, 377401575, 662824084, 1147902781, 2035015474]

lcg = LCG::Microsoft.new(1) p (1..5).map {lcg.rand}

1. prints [41, 18467, 6334, 26500, 19169]</lang>