home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / CPM / LANGUAGS / PASCAL / TAUSGEN3.PAS < prev    next >
Pascal/Delphi Source File  |  2000-06-30  |  12KB  |  392 lines
  1. program Tausgen ;
  2.  
  3. {
  4.  Demonstration of Tausworthe Pseudorandom Number Generator
  5.  30 April 1985
  6.  David P Kirschbaum, Toad Hall
  7.  7573 Jennings Lane, Fayetteville NC  28303
  8.  on ARPANet:  ABN.ISCAMS@USC-ISID
  9.  See extracts from Professor Thesen's paper (the source of this code) at
  10.  the end of this program.
  11.  
  12.  Expect an initial delay during the initialization of the byte table and
  13.  offset pointers when the program actually runs.  During actual implementation
  14.  in a program, the table would be initialized only once (until you wanted to
  15.  reseed again, anyway).
  16.  
  17.  If using my simple screen display version below, you have flow control.
  18.  Enter any key to pause, and enter RETURN to continue.  (Your system may
  19.  stop on its own occasionally if your system XOFFs high-speed printing to the
  20.  console.  Just enter RETURN to continue.)
  21.  
  22.  Enter Q or q to end.
  23. }
  24.  
  25. {Relax string parameter checking}
  26.  
  27. {V$-}
  28.  
  29. Type
  30.   Line  = string[79] ;
  31.   Menu  = string[25] ;
  32.  
  33.  
  34. const
  35.   pminusl = 97  ; {p = 98; q = 27}
  36.   maxbit  = 255 ;
  37.   maxs    = 26  ;
  38.   Selection : array[1..4] of Menu  =
  39.               ('1  Random Number Display ' ,
  40.                '2  Speed Run (no display)' ,
  41.                '3  Distribution Display  ' ,
  42.                '4  Quit                  ') ;
  43. var
  44.   Ch        : Char ;
  45.   Sentence  : Line ;
  46.   BlankLine : Line ;
  47.   Halted    : Boolean ;
  48.   Quit      : Boolean ;
  49.   Error     : Boolean ;
  50.   f,s       : Integer ;
  51.   x,y       : Integer ;
  52.   count     : Integer ;
  53.   prandom   : Integer ;
  54.   choice    : Integer ;
  55.   ByteTable : array[0..pminusl]of byte ;
  56.   RCounter  : array[0..maxbit] of Integer ;
  57.  
  58. procedure CenterPrint(Sentence : Line) ;
  59. var
  60.   Centered : Line ;
  61.  
  62.   begin
  63.     Centered := BlankLine ;
  64.     y := (80 - Length(Sentence)) div 2 ;
  65.     Insert(Sentence,Centered,y) ;
  66.     Writeln(Centered) ;
  67.   end ;
  68.  
  69.  
  70. procedure PlotXY(var number : integer; value : integer) ;
  71.  
  72.   begin
  73.     x := (number mod 12) * 4 + 12 ;
  74.     y := (number div 12) + 2 ;
  75.     GoToXY(x,y) ;
  76.     Write(value:4) ;
  77.   end ;  {PlotXY}
  78.  
  79.  
  80. procedure CheckUser ;
  81.  
  82.   begin
  83.     if KeyPressed then
  84.       begin
  85.         GoToXY(1,24) ;
  86.         Read(Ch) ;
  87.         Quit := Ch in ['Q', 'q'] ;
  88.       end ;  {KeyPressed}
  89.   end ; {CheckUser}
  90.  
  91. procedure ContPrompt ;
  92.  
  93.   begin
  94.     Writeln ;
  95.     GoToXY(1,22) ;
  96.     CenterPrint('Hit any key to continue:  ' ) ;
  97.     Writeln ;
  98.     GoToXY(54,22) ;
  99.     Repeat until KeyPressed = True ;
  100.   end ;
  101.  
  102. procedure TauInit(var f,s : integer) ;
  103. var
  104.   BitIndex  : integer ;
  105.   WorkTable : array[0..pminusl]of byte ;
  106.  
  107.   begin
  108.     for count := 0 to pminusl do
  109.       begin
  110.         ByteTable[count] := maxbit ;
  111.         WorkTable[count] := 0   ;
  112.       end ;  {count loop}
  113.     f := pminusl ;
  114.     s := maxs ;
  115.     Quit := False ;
  116.     for BitIndex := 1 to 16 do
  117.       begin
  118.         CheckUser ;
  119.         If NOT Quit then
  120.           begin
  121.             for count := 1 to 9800 do
  122.               begin
  123.                 if f < pminusl then
  124.                   f := f + 1
  125.                 else
  126.                   f := 0 ;
  127.                 if s < pminusl then
  128.                   s := s + 1
  129.                 else
  130.                   s := 0 ;
  131.                 ByteTable[f] := ByteTable[f] xor ByteTable[s] ;
  132.               end;  {count loop}
  133.             if BitIndex > 8 then
  134.               begin
  135.                 ClrScr ;
  136.                 Writeln ;
  137.                 Writeln('Loop ',BitIndex) ;
  138.                 for count := 0 to pminusl do
  139.                   begin
  140.                     if odd(ByteTable[count]) then
  141.                       WorkTable[count] := (WorkTable[count] div 2) + 128
  142.                     else
  143.                       WorkTable[count] := WorkTable[count] div 2 ;
  144.                     ByteTable[count] := ByteTable[count] div 2 ;
  145.                     PlotXY(count,WorkTable[count]) ;
  146.                   end;  {count loop}
  147.               end ; {BitIndex > 8}
  148.           end ;  {not Quit}
  149.       end ;  {BitIndex loop}
  150.     ClrScr ;
  151.     CenterPrint('Final Byte Table.') ;
  152.     for count := 0 to pminusl do
  153.       begin
  154.         ByteTable[count] := WorkTable[count] ;
  155.         PlotXY(count,ByteTable[count]) ;
  156.       end;  {count loop}
  157.     f := pminusl ;
  158.     s := maxs ;
  159.   end ;  {TauInit}
  160.  
  161.  
  162. function RByte : byte ;
  163.   begin
  164.     if f < pminusl then
  165.       f := f + 1
  166.     else
  167.       f := 0 ;
  168.     if s < pminusl then
  169.       s := s + 1
  170.     else
  171.       s := 0 ;
  172.     RByte := ByteTable[f] ;
  173.     ByteTable[f] := ByteTable[f] xor ByteTable[s] ;
  174.   end ;  {RByte}
  175.  
  176.  
  177. {main program               [Toad Hall]}
  178. begin
  179.   FillChar(BlankLine,79,chr(32)) ;
  180.   ClrScr ;
  181.   CenterPrint('A Toad Hall Demonstration.') ;
  182.   CenterPrint('of a Tausworthe Pseudorandom Number Generator.') ;
  183.   CenterPrint('as described by Thesen et al.') ;
  184.   Writeln ;
  185.   CenterPrint('Beginning table initialization.') ;
  186.   TauInit(f,s) ;
  187.   ContPrompt ;
  188.   ClrScr ;
  189.   Writeln ;
  190.   CenterPrint('Initialization complete.') ;
  191.   Writeln ;
  192.   for x := 1 to 4 do
  193.     CenterPrint(Selection[x]) ;
  194.   Writeln ;
  195.   Quit  := False ;
  196.   Error := True ;
  197.   repeat
  198.     GoToXY(30,10) ;
  199.     Write('Enter Selection (1-4):  ') ;
  200.     Readln(choice) ;
  201.     if choice in [1..4] then
  202.       Error := False ;
  203.     if choice = 4 then
  204.       Quit := True ;
  205.   until Error = False ;
  206.  
  207.   If NOT Quit then
  208.     begin
  209.       ClrScr ;
  210.       CenterPrint('Beginning random number generation.') ;
  211.       for prandom := 0 to maxbit do
  212.         RCounter[prandom] := 0 ;
  213.  
  214. { This is a screen display so you (human) can see the numbers.  Comment
  215.  this out, and use the next segment if you want clock timing.  Add in your
  216.  own file-related stuff if you want to collect lots of numbers for some sort
  217.  of statistical analysis.  Easy... but YOU can do it!
  218. }
  219.  
  220.       if choice = 1 then
  221.         repeat
  222.           for y := 5 to 23 do
  223.             begin
  224.               If NOT Quit then
  225.                 begin
  226.                   for x := 1 to 16 do
  227.                     begin
  228.                       GoToXY(x * 4,y) ;
  229.                       Write(RByte:4) ;
  230.                     end;  {x}
  231.                   CheckUser ;
  232.                 end ;  {NOT Quit}
  233.             end;  {y}
  234.         until Quit = True ;  {choice 1}
  235.  
  236. {And if you want to do some timing loops - this sucker is FAST!
  237.  Kludge in any system clock you have at the beginning and end of this
  238.  next routine, bring it on line instead of the display one above, and
  239.  do some timing.  On my Z80B CP/M system (Morrow Decision I), wristwatch
  240.  timing looked like 8 or 9 seconds for the 32000 loop.  Prof. Thesen
  241.  reported 8 seconds on an IBM-PC.}
  242.  
  243.         if choice = 2 then
  244.           begin
  245.             CenterPrint('30 000 random numbers.') ;
  246.             Writeln ;
  247.             CenterPrint('The big hand is pointing to .') ;
  248.             for count := 1 to 32000 do
  249.               prandom := RByte ;
  250.             CenterPrint('The big hand is pointing to .') ;
  251.             for x := 1 to 20 do
  252.               Write(' ') ;
  253.             Write('Completed with count = ',count) ;
  254.             Writeln(', and last number = ',prandom) ;
  255.           end ;  {choice 2}
  256.  
  257. {
  258.  And if you want a down-and-dirty analysis of random number distribution,
  259.  try this one ...
  260. }
  261.  
  262.         if choice = 3 then
  263.           repeat
  264.             CheckUser ;
  265.             prandom := RByte ;
  266.             if RCounter[prandom] < maxbit then
  267.               RCounter[prandom] := RCounter[prandom] + 1
  268.             else RCounter[prandom] := 0 ;
  269.             PlotXY(prandom,RCounter[prandom]) ;
  270.           until Quit = True ;  {choice 3}
  271.  
  272.     end  ;  {not Quit}
  273.  
  274.   GoToXY(1,23) ;
  275.   Writeln ;
  276.   Writeln ;
  277.   CenterPrint('Terminating.') ;
  278. end.
  279.  
  280. {
  281.  Algorithms extracted from
  282.  "Some Efficient Random Number Generators for Microcomputers"
  283.  by Arne Thesen, Shanshan Sun and Tzyh-Jong Wang
  284.  Industrial Engineering Department
  285.  University of Wisconsin-Madison
  286.  
  287.  (Released to the Public Domain by Professor Thesen, 23 April 1985
  288.  (message on file at Toad Hall))
  289.  
  290.  This is a Tausworthe Generator.  Extracts from Prof. Thesen's paper, the
  291.  implementing demonstration code, and some comments...
  292.  
  293.  by David Kirschbaum, Toad Hall (ABN.ISCAMS@USC-ISID
  294.  "A Tausworthe generator
  295.  
  296.  "Tausworthe [5] has suggested a different approach to the generation of
  297.  pseudo-random integers.  This procedure, which operates directly on bits
  298.  to form a stream of random bits, has been shown to produce random number
  299.  sequences that (1) have improved statistical properties over LC [linear
  300.  congruential] generators, and (2) have an arbitrarily long period independent
  301.  of the word size of the computer used.
  302.  
  303. "Tausworthe generators are not in widespread use.  This could be because they
  304. are difficult to implement efficiently in a higher level language, or because
  305. their improved statistical properties are of marginal utility on larger
  306. computers.  On microcomputers, the situation is quite different.  Here the
  307. improvement in period length and in statistical properties is quite
  308. substantial, and well written Tausworthe generators are not necessarily
  309. more time consuming than other classes of generators.
  310.  
  311. "Algorithm
  312.  
  313. The basic procedure of a Tausworthe type generator is illustrated in
  314. Figure 1:
  315.  
  316.          B[i-p]   B[i-r]     B[i]
  317.    B = . . x . . . . x . . . . x . . . . .
  318.            |         |         |
  319.            +-->xor<--+         |
  320.                 |              |
  321.                 +--------------+
  322.  
  323.  "Figure 1: Relationship between bits in a Tausworthe
  324.             generated bit stream.
  325.  
  326.  "Here B is defined as a sequence of bits, and the relationship between
  327.  individual bits in the sequence is defined in Algorithm 2:
  328.  
  329.     B[i] = B[i-r]  XOR  B[i-p]
  330.  
  331.     where:   i   = any integer
  332.          r , p   = fixed integers with 0<r<p.
  333.           XOR    = the exclusive OR operator
  334.                    yielding 0 if the terms are
  335.                    equal and 1 if they are not.
  336.  
  337.   "Algorithm 2: Tausworthe generator
  338.  
  339.  "When r and p are properly selected (as primitive trimodals [8]), the maximum
  340.  period of the stream (B) is 2**p - 1.
  341.  
  342.  "Implementation
  343.  
  344.  "In Program 3 we present a FORTRAN implementation of the Tausworthe based on
  345.  an idea first proposed by Lewis and Payne [3].  To avoid the need to access
  346.  individual bits, the algorithm maintains 15 independent and parallel streams
  347.  of bits, and the exclusive OR operation is performed on all bits at once.
  348.  Since each of the independent bit streams have a period of 2**p - 1, the
  349.  resulting stream of integers will also have a period of 2**p - 1."
  350.  
  351.  [Algorithm omitted - Toad Hall]
  352.  
  353.  "BIBLIOGRAPHY  [segment - Toad Hall]
  354.   ...
  355.   5.  Tausworthe, R.C., Random Numbers Generated by
  356.       Linear Recurrence, Modulo Two, Math. Comp. 19
  357.       (1965) 201-209.
  358.  
  359.   6.  Thesen, Arne, An Efficient Generator of Uniformly
  360.       Distributed Random Deviates Between Zero and One.
  361.       Technical Report.  Mathematics Research Center,
  362.       University of Wisconsin-Madison, 1983.  To appear
  363.       in Simulation.
  364.  
  365.   7.  Thesen, Arne and Tzyh-Jong Wang, Some Efficient
  366.       Random Number Generators for Micro Computers.  MRC
  367.       Technical Summary Report #2562.  Mathematics Research
  368.       Center, University of Wisconsin-Madison, 1983.
  369.  
  370.   8.  Zierler, Niel and John Brillhart, On Primitive
  371.       Trinomials (Mod 2), Information and Control,
  372.       Vol. 13 pp 541-554, 1968."
  373.  
  374.  Toad Hall Notes
  375.  
  376.  In the Pascal code from Prof. Thesen et al, a random byte (expressed as a
  377.  integer 0-255) is generated from eight parallel streams of random bits.  It
  378.  uses the xor operator built in to TurboPascal.
  379.  
  380.  The initial seed values pminusl and s were (I believe) chosen because of
  381.  their long period.  Not fully understanding all this yet, I expect there are
  382.  other suitable seeds (perhaps more of the primitive trimodals referred to in
  383.  the basic paper), or perhaps ANY value between 0 and 97 for pminusl, and any
  384.  lesser value for s?  (Sorry, guys, if my thinking is kind of dumb - I'm no
  385.  mathemetician - just a cut-and-try hacker.)
  386.  
  387.  Prof. Thesen states the Tausworthe algorithm is not constrained by the word
  388.  size of the computer.  This leads me to believe we might be able to implement
  389.  more than the 8 random bit streams he used in this example, but I expect it
  390.  will impact severely on the code speed.
  391. }
  392.