unit rmar;
{
 This random number generator originally appeared in "Toward a Universal
 Random Number Generator" by George Marsaglia and Arif Zaman.
 Florida State University Report: FSU-SCRI-87-50 (1987)

 It was later modified by F. James and published in "A Review of Pseudo-
 random Number Generators"

 THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
       (However, a newly discovered technique can yield
         a period of 10^600. But that is still in the development stage.)

 It passes ALL of the tests for random number generators and has a period
   of 2^144, is completely portable (gives bit identical results on all
   machines with at least 24-bit mantissas in the floating point
   representation).

 The algorithm is a combination of a Fibonacci sequence (with lags of 97
   and 33, and operation "subtraction plus one, modulo one") and an
   "arithmetic sequence" (using subtraction).
========================================================================
C language version was written by Jim Butler, and was based on a
FORTRAN program posted by David LaSalle of Florida State University.

Adapted for Delphi by Anton Zhuchkov (fireton@mail.ru) in February, 2002
}

interface


{ This is the initialization routine for the random number generator RANMAR()
  NOTE: The seed variables can have values between:    0 <= IJ <= 31328
                                                       0 <= KL <= 30081
  The random number sequences created by these two seeds are of sufficient
  length to complete an entire calculation with. For example, if sveral
  different groups are working on different parts of the same calculation,
  each group could be assigned its own IJ seed. This would leave each group
  with 30000 choices for the second seed. That is to say, this random
  number generator can create 900 million different subsequences -- with
  each subsequence having a length of approximately 10^30.

  Use IJ = 1802 & KL = 9373 to test the random number generator. The
  subroutine RANMAR should be used to generate 20000 random numbers.
  Then display the next six random numbers generated multiplied by 4096*4096
  If the random number generator is working properly, the random numbers
  should be:
            6533892.0  14220222.0  7275067.0
            6172232.0  8354498.0   10633180.0
}
procedure RMSeed(IJ, KL : Integer);

{  This is the random number generator proposed by George Marsaglia in
  Florida State University Report: FSU-SCRI-87-50
}
function RMRandom: Double;

const
	Seeded : Boolean = False;

implementation
uses SysUtils;

var
  u         : array [0..97] of Double;
	c, cd, cm : Double;
  i97, j97  : Integer;


procedure RMSeed(IJ, KL : Integer);
var
	i, j, k, l, ii, jj, m : Integer;
  s, t : Double;
begin
	if (ij<0) or (ij>31328) or (kl<0) or (kl>30081) then
  	raise Exception.Create('Random generator seed not within the valid range!');

	i := (ij div 177) mod 177 + 2;
	j := ij mod 177 + 2;
	k := (kl div 169) mod 178 + 1;
	l := kl mod 169;

	for ii:=1 to 97 do
  begin
		s := 0.0;
		t := 0.5;
		for jj:=1 to 24 do
    begin
			m := (((i*j) mod 179)*k) mod 179;
			i := j;
			j := k;
			k := m;
			l := (53*l + 1) mod 169;
			if ((l*m) mod 64 >= 32) then s := s + t;
			t := t*0.5;
		end;
		u[ii] := s;
	end;

	c := 362436.0 / 16777216.0;
	cd := 7654321.0 / 16777216.0;
	cm := 16777213.0 / 16777216.0;

	i97 := 97;
	j97 := 33;

	Seeded  := True;
end;



function RMRandom: Double;
var
	uni: Double;
begin
	if not Seeded then
  	raise Exception.Create('Random generator not seeded!');
  uni := u[i97] - u[j97];
  if (uni < 0.0) then uni := uni + 1.0;
  u[i97] := uni;
  Dec(i97);
  if i97 = 0 then i97 := 97;
  dec(j97);
  if j97 = 0 then j97 := 97;
  c := c - cd;
  if c < 0.0 then c := c + cm;
  uni := uni - c;
  if uni < 0.0 then uni := uni + 1.0;
  Result := uni;
end;

end.