A Random Number Generator using Preprocessor to emulate IAND, IEOR and REM

Copyright (c) Kednos Corporation 2008. All rights reserved
LAYERED PRODUCT:  Alpha PLI, V4.5          OP/SYS: VMS, V8.3
                  VAX   PLI, V3.8          OP/SYS  VMS, V7.3

SOURCE:     Kednos Customer Support Center


OVERVIEW:

Several statements have been commented out as they are not supported on 
OpenVMS.  The logical operators generate warning messages relating to the 
conversion from fixed bin the bit string.  Use of the BIT builtin will
remove these warnings.  They are left in here for instructive reasons.
Actually creating bit(32) strings BASED or DEFINED on these fixed bin
variables and applying the logical operators on these is more efficient.

/***********************************************************************/
/*                                                                     */
/* G. Marsaglia and W. W. Tsang,                                       */
/* "The 64-bit Universal RNG".                                         */
/* Generates double precision floating-point random numbers.           */
/*                                                                     */
/***********************************************************************/

/* Translated to PL/I from Fortran, by R. A. Vowels, 23 November 2005. */
/* modified to run on OpenVMS TLJ Linden 14-APR-2008                   */

/*(subscriptrange, fixedoverflow, size):*/
main: proc options (main);
   declare x float binary (53), i fixed binary (31);
   declare U(97) static float binary (53) external;

/* George Marsaglia & Wai Wan Tsang's 64-bit pseudo-random number generator */

   on fixedoverflow begin;
	end;
%IAND: proc(a,b) returns(char);
   dcl (a,b) char;
   return('bool('||a||','||b||','||'''0001'''||'B)');
   %end;

%IEOR: proc(a,b) returns(char);
   dcl (a,b) char;
   return('bool('||a||','||b||','||'''0110'''||'B)');
   %end;

%REM: proc(a,b) returns(char);
   dcl (a,b) char;
   return(a||'-'||b||'*trunc('||a||'/'||b||')');
   %end;

%ACTIVATE IEOR, IAND, REM;

uni64: procedure returns (float binary (53)) /*options(reorder)*/;
   declare r float binary (53) static /* initial
                        ((9007199254740881.0e0 / 9007199254740992.e0))*/;
   declare d float binary (53) static /*initial
                        ((362436069876.0e0 / 9007199254740992.0e0))*/;
   declare c float binary (53) static initial (0),
           i fixed binary  (7) static initial (97),
           j fixed binary  (7) static initial (33),
           x float binary (53);
   declare U(97) float binary (53) external;

   r = 9007199254740881.0e0 / 9007199254740992.e0;
   d = 362436069876.0e0 / 9007199254740992.0e0;

   x=U(i)-U(j); if x<0 then x=x+1.0; U(i)=x;
   i = i - 1; if i=0 then i=97; j = j - 1; if j=0 then j=97;
   c=c-d; if c<0 then c=c+r;
   x=x-c; if x<0 then return (x+1.); return (x); 
   end uni64;

/* A two-seed function for filling the static array U(97) one bit at a time. */
fillU: procedure (seed1, seed2) /* options(reorder)*/;
  declare (seed1, seed2) fixed binary (31);
  declare (s, t) float binary (53),
          (x, y) fixed binary (31),
	xb bit(32) based(addr(x)), yb bit(32) based(addr(y)),
          (i, j) fixed binary (7);
  declare U(97) float binary (53) external;

   x=seed1; y=seed2;
   do i=1 to 97;
      s=0.0; t=0.5;
      do j=1 to 53;
/*(nofixedoverflow, nosize):*/
         x=REM(6969*x, 65543);
/*(nofixedoverflow, nosize):*/
         y=REM(8888*y, 65579);
         if IAND(IEOR(xb,yb), 32 )>0 then s=s+t;
/*         if ((x ^ y) & bit(32,32)) > 0 then s=s+t; */
         t=.5*t;
      end;
      U(i)=s;
   end;
end fillU;

   call fillU (123456789,987654321);

   do i=1 to 10000000; x=uni64(); end;
   put skip list ('    x                 x*9007199254740992.0');
   do i=1 to 5;
      x = uni64();
      put skip edit(x, x*9007199254740992.0E0) (f(20,17), f(18));
   end;
/*
rem: proc(a,b)returns(fixed bin(31));
dcl (a,b) fixed bin(31);
return(a - b * trunc(a/b));
end rem;
*/
end main;