Linear congruential generator

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. $r_{n}$ and $r_{n+1}$ are not independent, as true random numbers would be. 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 Miller-Rabin primality test. 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.

In these formulas, the seed becomes $state_{0}$ . The random sequence is $rand_{1}$ , $rand_{2}$ and so on.

BSD formula:

$state_{n+1}=1103515245\times state_{n}+12345{\pmod {2^{31}}}$
$rand_{n}=state_{n}$

Microsoft formula:

$state_{n+1}=214013\times state_{n}+2531011{\pmod {2^{31}}}$
$rand_{n}=seed_{n}\div 2^{16}$

The BSD formula was so awful that FreeBSD switched to a different formula. More info is at Random number generator (included)#C.

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):

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; }

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

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

else /* MS rand */

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

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

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>

F#

This example may be incorrect.

These generators can yield negative numbers because the modulus operator of F# matches the sign of the first operand. The generators must only yield numbers from 0 to m - 1.

Negative numbers from lcg.bsd are congruent to the correct numbers: -740551042 is congruent to 1406932606. Negative numbers from lcg.ms are off by one, because the division truncated to the wrong direction: -11529 is congruent to 21239, but expected 21238.

Please verify it and remove this message. If the example does not match the requirements or does not work, replace this message with Template:incorrect or fix the code yourself.

<lang fsharp>module lcg =

   let bsd seed =
       let state = ref seed
       (fun (_:unit) ->
           state := (1103515245 * !state + 12345) % (1<<<31)
           !state)

   let ms seed =
       let state = ref seed
       (fun (_:unit) ->
           state := (214013 * !state + 2531011) % (1<<<31)
           !state / (1<<<16))

</lang>

let rnd = lcg.bsd 0; [for n in [0 .. 9] -> rnd()];;
BSD: [12345; -740551042; -1492899873; -698016724; 229283573; -1038148470; 1051550459; -853684456; -1353011855; 551188310]
MS: [38; 7719; -11529; 2437; -23912; 11797; 8365; 32285; -22317; 30612]

Fortran

Works with: Fortran version 90 and later

<lang fortran>module lcgs

 implicit none

 integer, parameter :: i64 = selected_int_kind(18)
 integer, parameter :: a1 = 1103515245, a2 = 214013
 integer, parameter :: c1 = 12345, c2 = 2531011
 integer, parameter :: div = 65536
 integer(i64), parameter :: m = 2147483648_i64  ! need to go to 64 bits because
                                                ! of the use of signed integers

contains

function bsdrand(seed)

 integer :: bsdrand
 integer, optional, intent(in) :: seed
 integer(i64) :: x = 0
 
 if(present(seed)) x = seed
 x = mod(a1 * x + c1, m)
 bsdrand = x

end function

function msrand(seed)

 integer :: msrand
 integer, optional, intent(in) :: seed
 integer(i64) :: x = 0

 if(present(seed)) x = seed 
 x = mod(a2 * x + c2, m)
 msrand = x / div

end function end module

program lcgtest

 use lcgs
 implicit none
 integer :: i
 
 write(*, "(a)") "      BSD            MS"
 do i = 1, 10
   write(*, "(2i12)") bsdrand(), msrand()
 end do

end program</lang> Output

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

Icon and Unicon

The following LCRNG's behave in the same way maintaining the state (seed) from round to round. The seed can be changed by specifying an argument. Otherwise it defaults to zero. <lang Icon>link printf

procedure main()

  printf("       BSD        MS\n")
  every 1 to 10 do 
     printf("%10s %10s\n",rand_BSD(),rand_MS())

end

procedure rand_BSD(x) static seed

  seed := \x | \seed | 0
  seed := (1103515245 * seed + 12345) % 2147483648
  return seed

end

procedure rand_MS(x) static seed

  seed := \x | \seed | 0
  seed := (214013 * seed + 2531011) % 2147483648
  return ishift(seed,-16)

end</lang>

Library: Icon Programming Library

printf.icn provides printf

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>

Perl 6

Works with: Niecza

Define subroutines implementing the LCG algorithm for each version then use those to generate lazy infinite lists of values and return the first 10 values from each. <lang perl6>my $mod = 2**31; sub bsd ($seed) { ( 1103515245 * $seed + 12345 ) % $mod }; sub ms ($seed) { ( 214013 * $seed + 2531011 ) % $mod };

say 'BSD LCG first 10 values:'; .say for ( 0.&bsd, -> $seed { $seed.&bsd } ... * )[^10];

say "\nMS LCG first 10 values:"; ($_ +> 16).say for ( 0.&ms, -> $seed { $seed.&ms } ... * )[^10];</lang>

BSD LCG first 10 values:
12345
1406932606
654583775
1449466924
229283573
1109335178
1051550459
1293799192
794471793
551188310

MS LCG first 10 values:
38
7719
21238
2437
8855
11797
8365
32285
10450
30612

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}

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

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

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

Tcl

Using an object-oriented solution, inspired by (but not a translation of) the Ruby solution above. <lang tcl>package require Tcl 8.6

General form of a linear-congruential RNG

oo::class create LCRNG {

   variable seed A B C D
   constructor {init a b c d} {

if {$init < 1} {set init [clock clicks]} variable seed $init A $a B $b C $c D $d

   }
   method rand {} {

set seed [expr {($A * $seed + $B) % $C}] return [expr {$seed / $D}]

   }
   method srand x {

set seed $x

}

Subclass to introduce constants

oo::class create BSDRNG {

   superclass LCRNG
   constructor Template:InitialSeed -1 {

next $initialSeed 1103515245 12345 [expr {2**31}] 1

} oo::class create MSRNG {

   superclass LCRNG
   constructor Template:InitialSeed -1 {

next $initialSeed 214013 2531011 [expr {2**31}] [expr {2**16}]

}</lang> Demo code: <lang tcl>proc sample rng {foreach - {1 2 3 4 5} {lappend r [$rng rand]}; join $r ", "} puts BSD:\t\[[sample [BSDRNG new 1]]\] puts MS:\t\[[sample [MSRNG new 1]]\]</lang> Output:

BSD:	[1103527590, 377401575, 662824084, 1147902781, 2035015474]
MS:	[41, 18467, 6334, 26500, 19169]