//==================================================== file = shuffle.c =====
//=  Program to shuffle a series of size N                                  =
//=   - Implements a paired shuffle (for M * N shuffles)                    =
//===========================================================================
//=  Notes:                                                                 =
//=    1) Input from input file "in.dat" to stdin (see example below)       =
//=        * Comments are bounded by "&" characters at the beginning and    =
//=          end of the comment block                                       =
//=    2) Outputs the series shuffled to stdout                             =
//=    3) Number of complete shuffles, M, is set in #define section         =
//=-------------------------------------------------------------------------=
//= Example "in.dat" file:                                                  =
//=                                                                         =
//=    & Sample series of data which can be integers or reals.              =
//=      There are 10 values in this file. &                                =
//=     1                                                                   =
//=     2                                                                   =
//=     3                                                                   =
//=     4                                                                   =
//=     5                                                                   =
//=     6                                                                   =
//=     7                                                                   =
//=     8                                                                   =
//=     9                                                                   =
//=    10                                                                   =
//=-------------------------------------------------------------------------=
//= Example output (for above "in.dat" and M = 5):                          =
//=                                                                         =
//=   7.000000                                                              =
//=   3.000000                                                              =
//=   5.000000                                                              =
//=   1.000000                                                              =
//=   4.000000                                                              =
//=   2.000000                                                              =
//=   9.000000                                                              =
//=   8.000000                                                              =
//=   10.000000                                                             =
//=   6.00000                                                               =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 shuffle.c                                                 =
//=-------------------------------------------------------------------------=
//=  Execute: shuffle < in.dat                                              =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (06/14/03) - Genesis (from shuffle1.c)                    =
//===========================================================================

//----- Include files -------------------------------------------------------
#include <stdio.h>                 // Needed for printf() and feof()
#include <stdlib.h>                // Needed for exit() and atof()
#include <string.h>                // Needed for strcmp()
#include <math.h>                  // Needed for ceil()

//----- Defines -------------------------------------------------------------
#define MAX_SIZE 2000000           // Maximum size of time series data array
#define M            100           // Block size for external shuffle

//----- Globals -------------------------------------------------------------
double     X[MAX_SIZE];            // Time series read from "in.dat"
int        N;                      // Number of values in "in.dat"

//----- Function prototypes -------------------------------------------------
void   load_X_array(void);         // Load X array
void   shuffle(void);              // Shuffle the X array
double rand_val(int seed);         // Jain's RNG

//===========================================================================
//=  Main program                                                           =
//===========================================================================
void main(void)
{
  int i;                           // Loop counter

  // Load the series X
  load_X_array();

  // Do the shuffle
  shuffle();

  // Output the shuffled series
  for (i=0; i<N; i++)
    printf("%f \n", X[i]);
}

//===========================================================================
//=  Function to load X array from stdin and determine N                    =
//===========================================================================
void load_X_array(void)
{
  char      temp_string[1024];     // Temporary string variable

  // Read all values into X
  N = 0;
  while(1)
  {
    scanf("%s", temp_string);
    if (feof(stdin)) goto end;

    // This handles a comment bounded by "&" symbols
    while (strcmp(temp_string, "&") == 0)
    {
      do
      {
        scanf("%s", temp_string);
        if (feof(stdin)) goto end;
      } while (strcmp(temp_string, "&") != 0);
      scanf("%s", temp_string);
      if (feof(stdin)) goto end;
    }

    // Enter value in array and increment array index
    X[N] = atof(temp_string);
    N++;

    // Check if MAX_SIZE data values exceeded
    if (N >= MAX_SIZE)
    {
      printf("*** ERROR - greater than %ld data values \n", MAX_SIZE);
      exit(1);
    }
  }

  // End-of-file escape
  end:

  return;
}

//===========================================================================
//=  Function to shuffle series X                                           =
//===========================================================================
void shuffle(void)
{
  int    from_index;              // From index
  int    to_index;                // To index
  double temp;                    // Temporary value
  int    i, j;                    // Loop counters

  // Seed the RNG
  rand_val(1);

  // Do a complete shuffle M times
  for (i=0; i<M; i++)
    for (j=0; j<N; j++)
    {
      // Determine random from_index and to_index
      from_index = j;
      to_index = (int) ceil(rand_val(0) * N) - 1.0;

      // Swap 'em
      temp = X[from_index];
      X[from_index] = X[to_index];
      X[to_index] = temp;
    }
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=========================================================================
double rand_val(int seed)
{
  const long  a =      16807;  // Multiplier
  const long  m = 2147483647;  // Modulus
  const long  q =     127773;  // m div a
  const long  r =       2836;  // m mod a
  static long x;               // Random int value
  long        x_div_q;         // x divided by q
  long        x_mod_q;         // x modulo q
  long        x_new;           // New x value

  // Set the seed if argument is non-zero and then return zero
  if (seed > 0)
  {
    x = seed;
    return(0.0);
  }

  // RNG using integer arithmetic
  x_div_q = x / q;
  x_mod_q = x % q;
  x_new = (a * x_mod_q) - (r * x_div_q);
  if (x_new > 0)
    x = x_new;
  else
    x = x_new + m;

  // Return a random value between 0.0 and 1.0
  return((double) x / m);
}