Linear congruential generator
The linear congruential generator is a very simple example of a random number generator. All linear congruential generators use this formula:
Where:
- is a seed.
- , , , ..., are the random numbers.
- , , are constants.
If one chooses the values of , and with care, then the generator produces a uniform distribution of integers from to .
LCG numbers have poor quality. and are not independent, as true random numbers would be. Anyone who knows can predict , 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 . 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 . The random sequence is , and so on.
BSD formula:
- is in range 0 to 2147483647.
Microsoft formula:
- is in range 0 to 32767.
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):
12345 1406932606 654583775 1449466924 229283573 1109335178 1051550459 1293799192 794471793 551188310 38 7719 21238 2437 8855 11797 8365 32285 10450 30612
bc
As with dc, bc has no bitwise operators. <lang bc>/* BSD rand */
define rand() { randseed = (randseed * 1103515245 + 12345) % 2147483648 return randseed }
randseed = 1 rand(); rand(); rand(); print "\n"
/* Microsoft rand */
define rand() { randseed = (randseed * 214013 + 2531011) % 2147483648 return randseed / 65536 }
randseed = 1 rand(); rand(); rand(); print "\n"</lang>
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>
dc
dc has no bitwise operations, so this program uses the modulus operator (2147483648 %
) and division (65536 /
). Fortunately, dc numbers cannot overflow to negative, so the modulus calculation involves only non-negative integers.
For BSD rand(): <lang dc>[*
* lrx -- (random number from 0 to 2147483647) * * Returns a number from the BSD rand() sequence. * Seeded by storing a seed in register R. *]sz
[lR 1103515245 * 12345 + 2147483648 % d sR]sr
[* Set seed to 1, then print the first 3 random numbers. *]sz 1 sR lrx psz lrx psz lrx psz</lang>
1103527590 377401575 662824084
For Microsoft rand(): <lang dc>[*
* lrx -- (random number from 0 to 32767) * * Returns a number from the Microsoft rand() sequence. * Seeded by storing a seed in register R. *]sz
[lR 214013 * 2531011 + 2147483648 % d sR 65536 /]sr
[* Set seed to 1, then print the first 3 random numbers. *]sz 1 sR lrx psz lrx psz lrx psz</lang>
41 18467 6334
F#
<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
<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>
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>
Liberty BASIC
<lang lb> 'by default these are 0 global BSDState global MSState
for i = 1 to 10
print randBSD()
next i
for i = 1 to 10
print randMS()
next i
function randBSD()
randBSD = (1103515245 * BSDState + 12345) mod (2 ^ 31) BSDState = randBSD
end function
function randMS()
MSState = (214013 * MSState + 2531011) mod (2 ^ 31) randMS = int(MSState / 2 ^ 16)
end function </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
Creates a magic scalar whose value is next in the LCG sequence when read.<lang perl>use strict; package LCG;
use overload '0+' => \&get;
use integer; sub gen_bsd { (1103515245 * shift() + 12345) % (1 << 31) }
sub gen_ms { my $s = (214013 * shift() + 2531011) % (1 << 31); $s, $s / (1 << 16) }
sub set { $_[0]->{seed} = $_[1] } # srand sub get { my $o = shift; ($o->{seed}, my $r) = $o->{meth}->($o->{seed}); $r //= $o->{seed} }
sub new { my $cls = shift; my %opts = @_; bless { seed => $opts{seed}, meth => $opts{meth} eq 'MS' ? \&gen_ms : \&gen_bsd, }, ref $cls || $cls; }
package main;
my $rand = LCG->new;
print "BSD:\n"; print "$rand\n" for 1 .. 10;
$rand = LCG->new(meth => 'MS');
print "\nMS:\n"; print "$rand\n" for 1 .. 10;</lang>output<lang>BSD: 12345 1406932606 654583775 1449466924 229283573 1109335178 1051550459 1293799192 794471793 551188310
MS: 38 7719 21238 2437 8855 11797 8365 32285 10450 30612</lang>
Perl 6
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]