program Tausgen ;

{
 Demonstration of Tausworthe Pseudorandom Number Generator
 30 April 1985
 David P Kirschbaum, Toad Hall
 7573 Jennings Lane, Fayetteville NC  28303
 on ARPANet:  ABN.ISCAMS@USC-ISID
 See extracts from Professor Thesen's paper (the source of this code) at
 the end of this program.

 Expect an initial delay during the initialization of the byte table and
 offset pointers when the program actually runs.  During actual implementation
 in a program, the table would be initialized only once (until you wanted to
 reseed again, anyway).

 If using my simple screen display version below, you have flow control.
 Enter any key to pause, and enter RETURN to continue.  (Your system may
 stop on its own occasionally if your system XOFFs high-speed printing to the
 console.  Just enter RETURN to continue.)

 Enter Q or q to end.
}

{Relax string parameter checking}

{V$-}

Type
  Line  = string[79] ;
  Menu  = string[25] ;


const
  pminusl = 97  ; {p = 98; q = 27}
  maxbit  = 255 ;
  maxs    = 26  ;
  Selection : array[1..4] of Menu  =
              ('1  Random Number Display ' ,
               '2  Speed Run (no display)' ,
               '3  Distribution Display  ' ,
               '4  Quit                  ') ;
var
  Ch        : Char ;
  Sentence  : Line ;
  BlankLine : Line ;
  Halted    : Boolean ;
  Quit      : Boolean ;
  Error     : Boolean ;
  f,s       : Integer ;
  x,y       : Integer ;
  count     : Integer ;
  prandom   : Integer ;
  choice    : Integer ;
  ByteTable : array[0..pminusl]of byte ;
  RCounter  : array[0..maxbit] of Integer ;

procedure CenterPrint(Sentence : Line) ;
var
  Centered : Line ;

  begin
    Centered := BlankLine ;
    y := (80 - Length(Sentence)) div 2 ;
    Insert(Sentence,Centered,y) ;
    Writeln(Centered) ;
  end ;


procedure PlotXY(var number : integer; value : integer) ;

  begin
    x := (number mod 12) * 4 + 12 ;
    y := (number div 12) + 2 ;
    GoToXY(x,y) ;
    Write(value:4) ;
  end ;  {PlotXY}


procedure CheckUser ;

  begin
    if KeyPressed then
      begin
        GoToXY(1,24) ;
        Read(Ch) ;
        Quit := Ch in ['Q', 'q'] ;
      end ;  {KeyPressed}
  end ; {CheckUser}

procedure ContPrompt ;

  begin
    Writeln ;
    GoToXY(1,22) ;
    CenterPrint('Hit any key to continue:  ' ) ;
    Writeln ;
    GoToXY(54,22) ;
    Repeat until KeyPressed = True ;
  end ;

procedure TauInit(var f,s : integer) ;
var
  BitIndex  : integer ;
  WorkTable : array[0..pminusl]of byte ;

  begin
    for count := 0 to pminusl do
      begin
        ByteTable[count] := maxbit ;
        WorkTable[count] := 0   ;
      end ;  {count loop}
    f := pminusl ;
    s := maxs ;
    Quit := False ;
    for BitIndex := 1 to 16 do
      begin
        CheckUser ;
        If NOT Quit then
          begin
            for count := 1 to 9800 do
              begin
                if f < pminusl then
                  f := f + 1
                else
                  f := 0 ;
                if s < pminusl then
                  s := s + 1
                else
                  s := 0 ;
                ByteTable[f] := ByteTable[f] xor ByteTable[s] ;
              end;  {count loop}
            if BitIndex > 8 then
              begin
                ClrScr ;
                Writeln ;
                Writeln('Loop ',BitIndex) ;
                for count := 0 to pminusl do
                  begin
                    if odd(ByteTable[count]) then
                      WorkTable[count] := (WorkTable[count] div 2) + 128
                    else
                      WorkTable[count] := WorkTable[count] div 2 ;
                    ByteTable[count] := ByteTable[count] div 2 ;
                    PlotXY(count,WorkTable[count]) ;
                  end;  {count loop}
              end ; {BitIndex > 8}
          end ;  {not Quit}
      end ;  {BitIndex loop}
    ClrScr ;
    CenterPrint('Final Byte Table.') ;
    for count := 0 to pminusl do
      begin
        ByteTable[count] := WorkTable[count] ;
        PlotXY(count,ByteTable[count]) ;
      end;  {count loop}
    f := pminusl ;
    s := maxs ;
  end ;  {TauInit}


function RByte : byte ;
  begin
    if f < pminusl then
      f := f + 1
    else
      f := 0 ;
    if s < pminusl then
      s := s + 1
    else
      s := 0 ;
    RByte := ByteTable[f] ;
    ByteTable[f] := ByteTable[f] xor ByteTable[s] ;
  end ;  {RByte}


{main program               [Toad Hall]}
begin
  FillChar(BlankLine,79,chr(32)) ;
  ClrScr ;
  CenterPrint('A Toad Hall Demonstration.') ;
  CenterPrint('of a Tausworthe Pseudorandom Number Generator.') ;
  CenterPrint('as described by Thesen et al.') ;
  Writeln ;
  CenterPrint('Beginning table initialization.') ;
  TauInit(f,s) ;
  ContPrompt ;
  ClrScr ;
  Writeln ;
  CenterPrint('Initialization complete.') ;
  Writeln ;
  for x := 1 to 4 do
    CenterPrint(Selection[x]) ;
  Writeln ;
  Quit  := False ;
  Error := True ;
  repeat
    GoToXY(30,10) ;
    Write('Enter Selection (1-4):  ') ;
    Readln(choice) ;
    if choice in [1..4] then
      Error := False ;
    if choice = 4 then
      Quit := True ;
  until Error = False ;

  If NOT Quit then
    begin
      ClrScr ;
      CenterPrint('Beginning random number generation.') ;
      for prandom := 0 to maxbit do
        RCounter[prandom] := 0 ;

{ This is a screen display so you (human) can see the numbers.  Comment
 this out, and use the next segment if you want clock timing.  Add in your
 own file-related stuff if you want to collect lots of numbers for some sort
 of statistical analysis.  Easy... but YOU can do it!
}

      if choice = 1 then
        repeat
          for y := 5 to 23 do
            begin
              If NOT Quit then
                begin
                  for x := 1 to 16 do
                    begin
                      GoToXY(x * 4,y) ;
                      Write(RByte:4) ;
                    end;  {x}
                  CheckUser ;
                end ;  {NOT Quit}
            end;  {y}
        until Quit = True ;  {choice 1}

{And if you want to do some timing loops - this sucker is FAST!
 Kludge in any system clock you have at the beginning and end of this
 next routine, bring it on line instead of the display one above, and
 do some timing.  On my Z80B CP/M system (Morrow Decision I), wristwatch
 timing looked like 8 or 9 seconds for the 32000 loop.  Prof. Thesen
 reported 8 seconds on an IBM-PC.}

        if choice = 2 then
          begin
            CenterPrint('30 000 random numbers.') ;
            Writeln ;
            CenterPrint('The big hand is pointing to .') ;
            for count := 1 to 32000 do
              prandom := RByte ;
            CenterPrint('The big hand is pointing to .') ;
            for x := 1 to 20 do
              Write(' ') ;
            Write('Completed with count = ',count) ;
            Writeln(', and last number = ',prandom) ;
          end ;  {choice 2}

{
 And if you want a down-and-dirty analysis of random number distribution,
 try this one ...
}

        if choice = 3 then
          repeat
            CheckUser ;
            prandom := RByte ;
            if RCounter[prandom] < maxbit then
              RCounter[prandom] := RCounter[prandom] + 1
            else RCounter[prandom] := 0 ;
            PlotXY(prandom,RCounter[prandom]) ;
          until Quit = True ;  {choice 3}

    end  ;  {not Quit}

  GoToXY(1,23) ;
  Writeln ;
  Writeln ;
  CenterPrint('Terminating.') ;
end.

{
 Algorithms extracted from
 "Some Efficient Random Number Generators for Microcomputers"
 by Arne Thesen, Shanshan Sun and Tzyh-Jong Wang
 Industrial Engineering Department
 University of Wisconsin-Madison

 (Released to the Public Domain by Professor Thesen, 23 April 1985
 (message on file at Toad Hall))

 This is a Tausworthe Generator.  Extracts from Prof. Thesen's paper, the
 implementing demonstration code, and some comments...

 by David Kirschbaum, Toad Hall (ABN.ISCAMS@USC-ISID
 "A Tausworthe generator

 "Tausworthe [5] has suggested a different approach to the generation of
 pseudo-random integers.  This procedure, which operates directly on bits
 to form a stream of random bits, has been shown to produce random number
 sequences that (1) have improved statistical properties over LC [linear
 congruential] generators, and (2) have an arbitrarily long period independent
 of the word size of the computer used.

"Tausworthe generators are not in widespread use.  This could be because they
are difficult to implement efficiently in a higher level language, or because
their improved statistical properties are of marginal utility on larger
computers.  On microcomputers, the situation is quite different.  Here the
improvement in period length and in statistical properties is quite
substantial, and well written Tausworthe generators are not necessarily
more time consuming than other classes of generators.

"Algorithm

The basic procedure of a Tausworthe type generator is illustrated in
Figure 1:

         B[i-p]   B[i-r]     B[i]
   B = . . x . . . . x . . . . x . . . . .
           |         |         |
           +-->xor<--+         |
                |              |
                +--------------+

 "Figure 1: Relationship between bits in a Tausworthe
            generated bit stream.

 "Here B is defined as a sequence of bits, and the relationship between
 individual bits in the sequence is defined in Algorithm 2:

    B[i] = B[i-r]  XOR  B[i-p]

    where:   i   = any integer
         r , p   = fixed integers with 0<r<p.
          XOR    = the exclusive OR operator
                   yielding 0 if the terms are
                   equal and 1 if they are not.

  "Algorithm 2: Tausworthe generator

 "When r and p are properly selected (as primitive trimodals [8]), the maximum
 period of the stream (B) is 2**p - 1.

 "Implementation

 "In Program 3 we present a FORTRAN implementation of the Tausworthe based on
 an idea first proposed by Lewis and Payne [3].  To avoid the need to access
 individual bits, the algorithm maintains 15 independent and parallel streams
 of bits, and the exclusive OR operation is performed on all bits at once.
 Since each of the independent bit streams have a period of 2**p - 1, the
 resulting stream of integers will also have a period of 2**p - 1."

 [Algorithm omitted - Toad Hall]

 "BIBLIOGRAPHY  [segment - Toad Hall]
  ...
  5.  Tausworthe, R.C., Random Numbers Generated by
      Linear Recurrence, Modulo Two, Math. Comp. 19
      (1965) 201-209.

  6.  Thesen, Arne, An Efficient Generator of Uniformly
      Distributed Random Deviates Between Zero and One.
      Technical Report.  Mathematics Research Center,
      University of Wisconsin-Madison, 1983.  To appear
      in Simulation.

  7.  Thesen, Arne and Tzyh-Jong Wang, Some Efficient
      Random Number Generators for Micro Computers.  MRC
      Technical Summary Report #2562.  Mathematics Research
      Center, University of Wisconsin-Madison, 1983.

  8.  Zierler, Niel and John Brillhart, On Primitive
      Trinomials (Mod 2), Information and Control,
      Vol. 13 pp 541-554, 1968."

 Toad Hall Notes

 In the Pascal code from Prof. Thesen et al, a random byte (expressed as a
 integer 0-255) is generated from eight parallel streams of random bits.  It
 uses the xor operator built in to TurboPascal.

 The initial seed values pminusl and s were (I believe) chosen because of
 their long period.  Not fully understanding all this yet, I expect there are
 other suitable seeds (perhaps more of the primitive trimodals referred to in
 the basic paper), or perhaps ANY value between 0 and 97 for pminusl, and any
 lesser value for s?  (Sorry, guys, if my thinking is kind of dumb - I'm no
 mathemetician - just a cut-and-try hacker.)

 Prof. Thesen states the Tausworthe algorithm is not constrained by the word
 size of the computer.  This leads me to believe we might be able to implement
 more than the 8 random bit streams he used in this example, but I expect it
 will impact severely on the code speed.
}