Back

unit uRand13;
{13 Dimensional Low Discrepancy Sobol Sequences
This unit realizes the  Lambda/Pi/Tau  sequence by V.M. Sobol.
The Numerical Monte-Carlo Methods, Moscow, Mir, 1979.
 The present version was written by E.B. Krissinel', (C) 1991.)
author: Eugene B. Krissinel http://www.fh.huji.ac.il/~krissinel/)
created: (February 02, 1991)
last modified: February 02, 2003)  }

interface

uses
  uMatTypes;

const
  Rand13Dim = 13;

type
  Rand13Point = array [1..Rand13Dim] of TFloat;
(*
  TRand13 = class
  private
    fQ: Rand13Point;
  public
    property Dim: IntType read GetDim write SetDim;
    property Q: Rand13Point read GetQ write SetQ;
  end;
*)
var
  Rand13Seed     : longint;
{ This is the current Random Generator's  state  variable;
  It may be saved and restored;
  The initial value should be assigned by  Randomize13. }

{   Call for the random initial point installation.}
procedure  Randomize13;

procedure Rand13(M: integer; var Q : Rand13Point );
{ M is the dimension of random-point space: 1 <= M <= Rand13Dim,
  Q is the current random point vector,
  it is placed in the first  M  elements of Q. }

{ ---------------------------------------------------------------------------- }

implementation

const
  Rand13Depth    = 20;
  NU : array [1..Rand13Depth,1..Rand13Dim] of longint =
  (
   (  524288,  524288,  524288,  524288,  524288,  524288,  524288,
      524288,  524288,  524288,  524288,  524288,  524288 ),
   (  262144,  786432,  262144,  786432,  262144,  786432,  262144,
      786432,  262144,  786432,  262144,  786432,  262144 ),
   (  131072,  655360,  917504,  131072,  655360,  917504,  393216,
      393216,  131072,  655360,  917504,  131072,  655360 ),
   (   65536,  983040,  720896,  327680,  196608,  458752,  720896,
      851968,  983040,   65536,  327680,  720896,  458752 ),
   (   32768,  557056,  425984, 1015808,  491520,  294912,  229376,
      753664,  360448,  950272,  163840,  819200,  622592 ),
   (   16384,  835584,  999424,  475136,  835584,  180224,  507904,
      933888,  475136,  999424,   49152,  245760,  835584 ),
   (    8192,  696320,  548864,  663552, 1024000,  647168,  106496,
       90112,  106496,  647168, 1040384,   73728,  417792 ),
   (    4096, 1044480,  323584,  602112,  577536,  258048,  659456,
      897024,  348160,  692224,  397312,  782336,  618496 ),
   (    2048,  526336,  952320,  886784,  362496,  395264,  632832,
      206848,  972800,  600064,   96256,  854016,  321536 ),
   (    1024,  789504,  738304,  152576,  777216,  723968,  922624,
      703488,   64512,  181248,  838656,  646144,  144384 ),
   (     512,  657920,  421376,  368128,  136704,  691712, 1027584,
      229888,  631296,  352768,  393680,  403968,  199168 ),
   (     256,  986880, 1047296,  945408,  470784,  378624,  335616,
      148224,   86272,  578304,  855296,  339712,  230656 ),
   (     128,  559232,  528000,  491648,  886912,  877696,  798080,
      483712,  841856,  281216,  432512,  369280,  753792 ),
   (      64,  838848,  265024,  737472, 1039424,   78528,  186560,
       45888, 1048512,  398400,  973888,  274880, 1032384 ),
   (      32,  699040,  919136,  532512,  530080, 1040864,  887840,
      568032,  795488,  731552,  688416,  928672,  385248 ),
   (      16, 1048560,  724976,  798800,  274224,  589808,  574160,
      327568,  397520, 1033424,  311408,  601648,  176880 ),
   (       8,   52296,  428040,  133368,  657656,  819208,  273800,
      189016,  882920,  123640,  123096,  500072,  723960 ),
   (       4,  786444, 1000452,  332916,  203916,  966668,  448708,
     1044588,  700052,  430564,  200948,  995596,  492908 ),
   (       2,  655370,  552462, 1031842,  505434,  516110,  756774,
      280906,  194774,  207722,  260870,  174986,   98978 ),
   (       1,  983055,  326411,  482707,  851901,  323591,  246491,
      658927,  453871,  247585,  790151,  716363,  574417 )
  );
  Base = 1048576.0;
  Rand13Limit = 1024*1024;     {  The same as Base  }

procedure Randomize13;
begin
  Randomize;
  Rand13Seed := round(Random*Rand13Limit);
end;

procedure Rand13(M : integer; var Q: Rand13Point);
var
  iQ     : array [1..Rand13Dim] of longint;
  i,j    : integer;
  N      : longint;

begin

  for i := 1 to M  do
    iQ[i] := 0;

  Rand13Seed := Rand13Seed mod Rand13Limit;
  N := Rand13Seed;

  i := 1;
  while N<>0  do
    begin
      if odd(N)  then
        for j:=1 to M  do
          iQ[j] := iQ[j] xor NU[i,j];
//      N := N div 2;
      N := N shr 1;
      inc(i);
    end;

  for i:=1 to M  do
    Q[i] := iQ[i]/Base;

  inc(Rand13Seed);

end;

end. {  of the  Rand13_   unit  }

Top

Back

Generated by Lore 's Source to HTML Converter ( http://www.newty.de/lsc/index.html  )