Advanced

Simple

Advanced

Simple

Algorithm creator(s)

Pierre L'Ecuyer

PB author(s)

Information Management Systems, Inc.

Description

An excellent composite pseudo-RNG with an extremely long period ( >8*10^12 )

Note

Source: Pierre L'Ecuyer: Communications of the ACM, 31, 6, June 88

Source

http://www.infoms.com/bsp/list7-10.zip

See also

• Information Management Systems, Inc.

Source Code

Download source code file lecuyer.bas (Right-click -> "Save as ...")

DECLARE SUB SeedStore(Operation$, seed1%, seed2%, seed3%)

' RanSeed
' Seed the L'Ecuyer LCG3 Pseudo Random Number Generator

' Input: Seed1%  Seed value generator 1. Value: 1..32362 (inclusive)
'        Seed2%  Seed value generator 2. Value: 1..31726 (inclusive)
'        Seed3%  Seed value generator 3. Value: 1..31656 (inclusive)

SUB RanSeed(Seed1%, Seed2%, Seed3%) Public

  '*** Force lowerlimits to prevent lockup
  IF Seed1% < 1 THEN Seed1%=1
  IF Seed2% < 1 THEN Seed2%=1
  IF Seed3% < 1 THEN Seed3%=1

  '*** Force upperlimits to prevent wrong max. seeds
  IF Seed1%>32362 THEN Seed1%=32362
  IF Seed2%>31726 THEN Seed2%=31726
  IF Seed3%>31656 THEN Seed3%=31656

  '*** Store seeds
  SeedStore "S", Seed1%, Seed2%, Seed3%

END SUB 'RanSeed

' RanGetSeed
' Read the seedvalues. Read-Only

SUB RanGetSeed(Seed1%, Seed2%, Seed3%) Public

  SeedStore "R", Seed1%, Seed2%, Seed3%

END SUB 'RanGetSeed

' Rand
' Composite Pseudo Random Number Generator
' Source: Pierre L'Ecuyer, Communications of the ACM, 31, 6, June 88
' An excellent Pseudo RNG with an extremely long period ( >8*10^12).

' Before using this routine, initialize by calling RanSeed

' Returns: single precision uniformly distributed random number
' from range  0.0 ... 1.0

FUNCTION Rand!() Public

  LOCAL seed1%, seed2%, seed3%, itemp%, icombined%

  ' *** Read current seeds
  SeedStore "R", seed1%, seed2%, seed3%

  '*** First algorithm
  itemp% = seed1% \ 206
  seed1% = 157 * (seed1% - itemp%*206) - itemp%*21
  IF seed1% < 0 THEN seed1% = seed1% + 32363

  '*** Second algorithm
  itemp% = seed2% \ 217
  seed2% = 146 * (seed2% - itemp%*217) -itemp%*45
  IF seed2% < 0 THEN seed2% = seed2% + 31727

  '*** Third algorithm
  itemp%= seed3% \ 222
  seed3%= 142 * (seed3% - itemp%*222) - itemp%*133
  IF seed3% <0 THEN seed3% = seed3% + 31657

  '*** Combine the random number generators
  icombined% = seed1% - seed2%
  IF icombined% > 706 THEN icombined% = icombined% - 32362
  icombined% = icombined% + seed3%
  IF icombined% < 1 THEN icombined% = icombined% + 32362

  '*** Scale to single precision value (0.0 , 1.0]
  Rand! = icombined% * 3.0899E-5

  '*** Update internal seeds
  SeedStore "S", seed1%, seed2%, seed3%

END FUNCTION 'Rand

' SeedStore
' Store and Retrieve seeds
'   Input: Operation$ "R"=Retrieve seedvalues, "S"=Store seedvalues
'   Input/output: Seed1%, Seed2%, Seed3%

SUB SeedStore(Operation$, seed1%, seed2%, seed3%) Private
  STATIC static_seed1%, static_seed2%, static_seed3%, static_init?

  SELECT CASE UCASE$(Operation$)
    CASE "S"              'Store seed values
      static_init?  = 1      'Set initialization flag
      static_seed1% = seed1%
      static_seed2% = seed2%
      static_seed3% = seed3%
    CASE "R"              'Retrieve seed values
      IF static_init? = 1 THEN   'Initialized; use seeds
        seed1% = static_seed1%
        seed2% = static_seed2%
        seed3% = static_seed3%
      ELSE                   ' Not initialized. Force random seeds
        RANDOMIZE TIMER
        static_seed1% = 1 + INT( RND * 32362 )
        static_seed2% = 1 + INT( RND * 31726 )
        static_seed3% = 1 + INT( RND * 31656 )
        static_init?  = 1       ' Consider now as initialized

        seed1% = static_seed1%
        seed2% = static_seed2%
        seed3% = static_seed3%
      END IF
  END SELECT

END SUB 'SeedStore
 
DECLARE SUB RanPermutation(Arr%(), NrElements%)

'FUNCTION PBMAIN  ' uncomment for use with PB/CC

DIM Arr%(1:10)
CLS

' *** Initialize array of 10 elements
FOR i%=1 TO 10
  Arr%(i%)=i%
NEXT i

' *** Initialize the L'Ecuyer LCG3 random number generator
Seed1%=12345
Seed2%=12345
Seed3%=12345
RanSeed Seed1%, Seed2%, Seed3%

' *** Use the L'Ecuyer RNG to produce a random permutation
RanPermutation Arr%(), 10

' *** Print the random permutation
PRINT "Random permutation: ";
FOR i%=1 TO 10
  PRINT Arr%(i%) " ";
NEXT i%

' *** Print the current internal seeds
RanGetseed seed1%, seed2%, seed3%
PRINT
PRINT "Current seeds: "seed1%, seed2%, seed3%

END 'FUNCTION  ' uncomment for use with PB/CC


SUB RanPermutation(Arr%(), NrElements%) PUBLIC
  LOCAL i%, swapelement%

  FOR i% = 1 TO NrElements%-1
    swapelement% = INT(Rand! * (NrElements%-i%+1))
    SWAP Arr%(i%), Arr%(i%+swapelement%)
  NEXT i

END SUB 'RanPermutation