(50G) Quasi Monte Carlo Halton Sequence Programme

[Gerald H]

Edit 2017-01-28: Names of variables in documentation updated to conform to programme.

This group of programmes generates a Quasi-
Monte Carlo series of D-dimension
Halton sequence vectors.

The problem with (quasi) random number generators is they clump: Monte Carlo sequences don't.

QMC is useful in, eg

https://en.wikipedia.org/wiki/Quasi-Monte_Carlo_method

Halton sequence ref:

https://en.wikipedia.org/wiki/Halton_sequence

DINIT stores the vector dimension
in D,

NINIT a precision parameter in
N &

‡INIT a list of primes in ‡LIST.

D must be =< length of ‡LIST.

N is the maximum ordinality
of the series.

QMCSeed takes an integer input
fixing the index in the series &
sets Qij & Dij, two lists of
lists of parameters, & x,the
current Halton sequence vector.

QMC takes no input & returns the
next vector in the series, also
stored in x.

Code:


QMC

::
  CK0
  ID ‡LIST
  CKLIST_
  ID x
  CKLIST_
  ID Qij
  CKLIST_
  ID Dij
  CKLIST_
  ID D
  CKREAL
  COERCEDUP
  ID N
  CKREAL
  %1+
  %LN
  BINT0
  FPTR2 ^3LAMBIND
  #1+_ONE_DO
  OVERINDEX@
  NTHCOMPDROP
  CKLIST_
  OVERINDEX@
  NTHCOMPDROP
  CKLIST_
  5PICK
  INDEX@
  NTHCOMPDROP
  CKREAL
  2GETLAM
  8PICK
  INDEX@
  NTHCOMPDROP
  CKREAL
  DUP1PUTLAM_
  %LN
  %/
  %CEIL
  COERCE
  #2+
  BINT2
  DO
  3PICK
  INDEX@
  NTHCOMPDROP
  CKREAL
  %+SWAP
  INDEX@
  #1-
  2DUP
  NTHCOMPDROP
  CKREAL
  %1+
  DUP4UNROLL
  SWAPROT
  PUTLIST
  3UNROLL
  1GETLAM
  %<
  ITE
  ExitAtLOOP
  ::
    %0
    INDEX@
    #1-
    4ROLL
    PUTLIST
    SWAP3PICK
    INDEX@
    #1-
    NTHCOMPDROP
    CKREAL
    %-
  ;
  LOOP
  ROT
  INDEX@
  6ROLL
  PUTLIST
  ROT
  INDEX@
  5ROLL
  PUTLIST
  4ROLL
  3UNROLL
  5ROLL
  4UNROLL
  LOOP
  '
  ID Dij
  ?STO_HERE
  '
  ID Qij
  ?STO_HERE
  2DROP
  3GETLAM
  ABND
  {}N
  DUP'
  ID x
  ?STO_HERE
;

QMCSeed 

::
  CK1&Dispatch
  BINT1
  ::
    ID ‡LIST
    CKLIST_
    ID D
    CKREAL
    COERCE
    ID N
    CKREAL
    %1+
    %LN
    OVER
    FPTR2 ^2LAMBIND
    NULL{}
    NULL{}
    5ROLL
    4ROLL
    #1+_ONE_DO
    {
      %1
    }
    NULL{}
    ROTDUP
    7PICK
    INDEX@
    NTHCOMPDROP
    CKREAL
    %0
    %1
    2GETLAM
    4PICK
    %LN
    %/
    %CEIL
    COERCE
    ZERO_DO
    3PICK
    %/
    7ROLL
    OVER
    >TCOMP
    7UNROLL
    4PICK
    4PICK
    %MOD
    7ROLL
    OVER
    >TCOMP
    7UNROLL
    5ROLL
    OVER
    %-
    5PICK
    %/
    5UNROLL
    OVER
    %*
    ROT
    %+
    SWAPLOOP
    SWAP
    5UNROLL
    3DROP
    5ROLL
    4ROLL
    >TCOMP
    5UNROLL
    4ROLL
    4ROLL
    >TCOMP
    4UNROLL
    SWAP
    5UNROLL
    LOOP
    DROP
    '
    ID Dij
    ?STO_HERE
    '
    ID Qij
    ?STO_HERE
    DROP
    1GETABND
    {}N
    '
    ID x
    ?STO_HERE
  ;
;

NINIT 

::
  CK1
  FPTR2 ^CK1Z
  FPTR2 ^Z>R
  '
  ID N
  ?STO_HERE
;

DINIT

::
  CK1
  FPTR2 ^CK1Z
  FPTR2 ^Z>R
  '
  ID D
  ?STO_HERE
;

‡INIT 

::
  CK0
  {
    %2
    %3
    %5
    %7
    %11
    %13
    %17
    %19
    %23
    % 29.
  }
  '
  ID ‡LIST
  ?STO_HERE
;