/***************************************************/
/* Program: shuffle.c                              */
/* Date:    091425PST Jan 1997                     */
/* Name:    Alfred A. Aburto Jr.                   */
/* Email:   aburto@nosc.mil, aburto@cts.com        */
/* Place:   San Diego, CA, USA                     */
/*                                                 */
/* This very simple routine tests a random number  */
/* generator using a card shuffling procedure.     */
/*                                                 */
/* Shuffle.c counts the number of occurrences of   */
/* card type (ace, two, three, ..., queen, king)   */
/* at position 103 in a stack of cards ( 4 decks   */
/* of 52 cards) for 26,000 shuffles. Ideally the   */
/* number of cards observed of the same type at    */
/* position 103 (or any other position) would be   */
/* 26000 / 13 = 2000 assuming the random number    */
/* generator provides a uniform distribution.      */
/***************************************************/

/***************************************************************/
/* Timer options. You MUST uncomment one of the options below  */
/* or compile, for example, with the '-DUNIX' option.          */
/***************************************************************/
/* #define Amiga       */
/* #define UNIX        */
/* #define UNIX_Old    */
/* #define VMS         */
/* #define BORLAND_C   */
/* #define MSC         */
/* #define MAC         */
/* #define IPSC        */
/* #define FORTRAN_SEC */
/* #define GTODay      */
/* #define CTimer      */
/* #define UXPM        */
/* #define MAC_TMgr    */
/* #define PARIX       */
/* #define POSIX       */
/* #define WIN32       */
/* #define POSIX1      */
/***********************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/****************************************/
/* Parameters for the mother of random  */
/* number generation programs: Mother() */
/****************************************/
#define m16Long   65536L        /* 2^16 */
#define m16Mask   0xFFFF        /* mask for lower 16 bits */
#define m15Mask   0x7FFF        /* mask for lower 15 bits */
#define m31Mask   0x7FFFFFFF    /* mask for 31 bits */
#define m32Double 4294967295.0  /* 2^32-1 */

static short mother1[10];
static short mother2[10];
static short mStart=1;

#define NDecks    4
int card[52*NDecks];
int suit[14];

main()
{
 int  i, j, k, l, m, n;   
 int  ICard, NCards, NShuf, shuffle();
 
 double x, y, z;
 double starttime, benchtime, dtime();

 char cname[13][10]  = { "ACE", "TWO", "THREE",
			 "FOUR", "FIVE", "SIX",
			 "SEVEN", "EIGHT", "NINE",
			 "TEN", "JACK", "QUEEN",
			 "KING"
			 };

 /*####################*/
 /*#  Initialization  #*/
 /*####################*/
 NCards = 52 * NDecks;
 NShuf  = 26000;
 ICard  = 103;

 for (i = 0; i < 14; i++)
    {
     suit[i] = 0;
    }

 printf("SHUFFLE: Test of Random Number Generator.\n");
 printf("Counts number of occurrences of card type at position\n");
 printf("%d in card deck for %d shuffles. Ideal would be\n",ICard,NShuf);
 printf("%d occurrences for each type of card.\n", NShuf/13);
 printf(" Card   Number    Percent\n");
 printf(" value  observed  from ideal\n");
 
 /********************************/
 /* Preset the NDecks of cards.  */
 /*  1 is Ace                    */
 /* 11 is Jack                   */
 /* 12 is Queen                  */
 /* 13 is King                   */
 /* We ignore the suit...        */
 /********************************/
 starttime = dtime();

 for (k = 0; k < NDecks; k++)
    {
     n = 52 * k;
     for (j = 0; j < 4; j++)
	  {
	   m = 13 * j + n;
	   for (i = 1; i < 14; i++)
		{
		 card[i + m - 1] = i;
		}
	  }
    }

 /***********************************/
 /* Pre-shuffle the NDecks of cards */
 /***********************************/
 shuffle();

 for (i = 1; i <= NShuf; i++)
    {
     shuffle();
     j = card[ICard];
     suit[j] = suit[j] + 1;
    }
 
 j = 0;
 y = 13.0 / (double)NShuf;
 for (i = 1; i < 14; i++)
    {
     k = suit[i];
     x = 100.0 * ((double)k * y - 1.0);
     printf("%5s   %5d    %5.1f\n",cname[i-1],k,x);
     j = j + k;
    }

 benchtime = dtime() - starttime;
 
 printf("Total Cards = %d\n",j);
 printf("Elapsed Time (sec): %f\n",benchtime);

 exit(0);
}

int shuffle()
{
 int i, j, k, l, m, n;
 unsigned long kr = 3141592;
 double r, x, scale, Mother();
 
 /**************************************/
 /*  If you want to test rand() you    */
 /*  will need to set the value of     */
 /*  scale as follows:                 */
 /*  scale  = 1.0 / (double)RAND_MAX;  */
 /*                                    */
 /*  For Mother() however scale = 1.0; */
 /**************************************/
 scale = 1.0;

 j = 52 * NDecks;
 r = (double)j * scale;
 k = 10 * j;
 for (i = 1; i <= k; i++)
    {
     m = (int)(r * Mother(&kr));
     n = (int)(r * Mother(&kr));
     /****************************/
     /*  m = (int)(r * rand());  */
     /*  n = (int)(r * rand());  */
     /****************************/
     l = card[n];
     card[n] = card[m];
     card[m] = l;
    }

 return 0;
}

/* Mother *************************************************************
|  George Marsaglia's The mother of all random number generators
|   producing uniformly distributed pseudo random 32 bit values
|   with period about 2^250.
|  The text of Marsaglia's posting is provided in the file mother.txt
|
|  The arrays mother1 and mother2 store carry values in their
|   first element, and random 16 bit numbers in elements 1 to 8.
|   These random numbers are moved to elements 2 to 9 and a new
|   carry and number are generated and placed in elements 0 and 1.
|  The arrays mother1 and mother2 are filled with random 16 bit
|   values on first call of Mother by another generator.  mStart
|   is the switch.
|
|  Returns:
|   A 32 bit random number is obtained by combining the output of the
|   two generators and returned in *pSeed.  It is also scaled by
|   2^32-1 and returned as a double between 0 and 1
|
|  SEED:
|   The inital value of *pSeed may be any long value
|
|   Bob Wheeler 8/8/94
*/

double Mother(unsigned long *pSeed)
{
 unsigned long  number, number1, number2;
 short n, *p;
 unsigned short sNumber;

 /* Initialize motheri with 9 random values the first time */
 if (mStart) 
   {
    sNumber= *pSeed&m16Mask;   /* The low 16 bits */
    number = *pSeed&m31Mask;   /* Only want 31 bits */

    p = mother1;
    for (n=18;n--;)
	 {
	  number=30903*sNumber+(number>>16);   /* One line multiply-with-cary */
	  *p++=sNumber=number&m16Mask;
	  if (n==9) p = mother2;
	 }
    
    /* make cary 15 bits */
    mother1[0]&=m15Mask;
    mother2[0]&=m15Mask;
    mStart=0;
   }

 /* Move elements 1 to 8 to 2 to 9 */
 memmove(mother1+2,mother1+1,8*sizeof(short));
 memmove(mother2+2,mother2+1,8*sizeof(short));

 /* Put the carry values in numberi */
 number1=mother1[0];
 number2=mother2[0];

 /* Form the linear combinations */
 number1+=1941*mother1[2]+1860*mother1[3]+1812*mother1[4]+
 1776*mother1[5]+1492*mother1[6]+1215*mother1[7]+
 1066*mother1[8]+12013*mother1[9];
	
 number2+=1111*mother2[2]+2222*mother2[3]+3333*mother2[4]+
 4444*mother2[5]+5555*mother2[6]+6666*mother2[7]+
 7777*mother2[8]+9272*mother2[9];

 /* Save the high bits of numberi as the new carry */
 mother1[0]=number1/m16Long;
 mother2[0]=number2/m16Long;
		
 /* Put the low bits of numberi into motheri[1] */
 mother1[1]=m16Mask&number1;
 mother2[1]=m16Mask&number2;

 /* Combine the two 16 bit random numbers into one 32 bit */
 *pSeed=(((long)mother1[1])<<16)+(long)mother2[1];

 /* Return a double value between 0 and 1 */
 return ((double)*pSeed)/m32Double;
}

/*****************************************************/
/* Various timer routines.                           */
/* Al Aburto, aburto@nosc.mil, 18 Feb 1997           */
/*                                                   */
/* t = dtime() outputs the current time in seconds.  */
/* Use CAUTION as some of these routines will mess   */
/* up when timing across the hour mark!!!            */
/*                                                   */
/* For timing I use the 'user' time whenever         */
/* possible. Using 'user+sys' time is a separate     */
/* issue.                                            */
/*                                                   */
/* Example Usage:                                    */
/* [timer options added here]                        */
/* main()                                            */
/* {                                                 */
/* double starttime,benchtime,dtime();               */
/*                                                   */
/* starttime = dtime();                              */ 
/* [routine to time]                                 */
/* benchtime = dtime() - starttime;                  */
/* }                                                 */
/*                                                   */
/* [timer code below added here]                     */
/*****************************************************/

/*********************************/
/* Timer code.                   */
/*********************************/
/*******************/
/*  Amiga dtime()  */
/*******************/
#ifdef Amiga
#include <ctype.h>
#define HZ 50

double dtime()
{
 double q;

 struct tt
       {
	long  days;
	long  minutes;
	long  ticks;
       } tt;

 DateStamp(&tt);

 q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;

 return q;
}
#endif

/*****************************************************/
/*  UNIX dtime(). This is the preferred UNIX timer.  */
/*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
/*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
/*****************************************************/
#ifdef UNIX
#include <sys/time.h>
#include <sys/resource.h>

#ifdef hpux
#include <sys/syscall.h>
#define getrusage(a,b) syscall(SYS_getrusage,a,b)
#endif

struct rusage rusage;

double dtime()
{
 double q;

 getrusage(RUSAGE_SELF,&rusage);

 q = (double)(rusage.ru_utime.tv_sec);
 q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
	
 return q;
}
#endif

/***************************************************/
/*  UNIX_Old dtime(). This is the old UNIX timer.  */
/*  Make sure HZ is properly defined in param.h !! */
/***************************************************/
#ifdef UNIX_Old
#include <sys/types.h>
#include <sys/times.h>
#include <sys/param.h>

#ifndef HZ
#define HZ 60
#endif

struct tms tms;

double dtime()
{
 double q;

 times(&tms);

 q = (double)(tms.tms_utime) / (double)HZ;
	
 return q;
}
#endif

/*********************************************************/
/*  VMS dtime() for VMS systems.                         */
/*  Provided by: RAMO@uvphys.phys.UVic.CA                */
/*  Some people have run into problems with this timer.  */
/*********************************************************/
#ifdef VMS
#include time

#ifndef HZ
#define HZ 100
#endif

struct tbuffer_t
      {
       int proc_user_time;
       int proc_system_time;
       int child_user_time;
       int child_system_time;
      };

struct tbuffer_t tms;

double dtime()
{
 double q;

 times(&tms);

 q = (double)(tms.proc_user_time) / (double)HZ;
	
 return q;
}
#endif

/******************************/
/*  BORLAND C dtime() for DOS */
/******************************/
#ifdef BORLAND_C
#include <ctype.h>
#include <dos.h>
#include <time.h>

#define HZ 100
struct time tnow;

double dtime()
{
 double q;

 gettime(&tnow);

 q = 60.0 * (double)(tnow.ti_min);
 q = q + (double)(tnow.ti_sec);
 q = q + (double)(tnow.ti_hund)/(double)HZ;
	
 return q;
}
#endif

/**************************************/
/*  Microsoft C (MSC) dtime() for DOS */
/**************************************/
#ifdef MSC
#include <time.h>
#include <ctype.h>

#define HZ CLOCKS_PER_SEC
clock_t tnow;

double dtime()
{
 double q;

 tnow = clock();

 q = (double)tnow / (double)HZ;
	
 return q;
}
#endif

/*************************************/
/*  Macintosh (MAC) Think C dtime()  */
/*************************************/
#ifdef MAC
#include <time.h>

#define HZ 60

double dtime()
{
 double q;

 q = (double)clock() / (double)HZ;
	
 return q;
}
#endif

/************************************************************/
/*  iPSC/860 (IPSC) dtime() for i860.                       */
/*  Provided by: Dan Yergeau, yergeau@gloworm.Stanford.EDU  */
/************************************************************/
#ifdef IPSC
extern double dclock();

double dtime()
{
 double q;

 q = dclock();
	
 return q;
}
#endif

/**************************************************/
/*  FORTRAN dtime() for Cray type systems.        */
/*  This is the preferred timer for Cray systems. */
/**************************************************/
#ifdef FORTRAN_SEC

fortran double second();

double dtime()
{
 double q;

 second(&q);
	
 return q;
}
#endif

/***********************************************************/
/*  UNICOS C dtime() for Cray UNICOS systems.  Don't use   */
/*  unless absolutely necessary as returned time includes  */
/*  'user+system' time.  Provided by: R. Mike Dority,      */
/*  dority@craysea.cray.com                                */
/***********************************************************/
#ifdef CTimer
#include <time.h>

double dtime()
{
 double    q;
 clock_t   clock(void);

 q = (double)clock() / (double)CLOCKS_PER_SEC;

 return q;
}
#endif

/********************************************/
/* Another UNIX timer using gettimeofday(). */
/* However, getrusage() is preferred.       */
/********************************************/
#ifdef GTODay
#include <sys/time.h>

struct timeval tnow;

double dtime()
{
 double q;

 gettimeofday(&tnow,NULL);
 q = (double)tnow.tv_sec + (double)tnow.tv_usec * 1.0e-6;

 return q;
}
#endif

/*****************************************************/
/*  Fujitsu UXP/M timer.                             */
/*  Provided by: Mathew Lim, ANUSF, M.Lim@anu.edu.au */
/*****************************************************/
#ifdef UXPM
#include <sys/types.h>
#include <sys/timesu.h>
struct tmsu rusage;

double dtime()
{
 double q;

 timesu(&rusage);

 q = (double)(rusage.tms_utime) * 1.0e-06;
	
 return q;
}
#endif

/**********************************************/
/*    Macintosh (MAC_TMgr) Think C dtime()    */
/*   requires Think C Language Extensions or  */
/*    #include <MacHeaders> in the prefix     */
/*  provided by Francis H Schiffer 3rd (fhs)  */
/*         skipschiffer@genie.geis.com        */
/**********************************************/
#ifdef MAC_TMgr
#include <Timer.h>
#include <stdlib.h>

static TMTask   mgrTimer;
static Boolean  mgrInited = false;
static double   mgrClock;

#define RMV_TIMER RmvTime( (QElemPtr)&mgrTimer )
#define MAX_TIME  1800000000L
/* MAX_TIME limits time between calls to */
/* dtime( ) to no more than 30 minutes   */
/* this limitation could be removed by   */
/* creating a completion routine to sum  */
/* 30 minute segments (fhs 1994 feb 9)   */

static void Remove_timer( )
{
 RMV_TIMER;
 mgrInited = false;
}

double  dtime( )
{
 if( mgrInited ) {
   RMV_TIMER;
   mgrClock += (MAX_TIME + mgrTimer.tmCount)*1.0e-6;
 } else {
   if( _atexit( &Remove_timer ) == 0 ) mgrInited = true;
   mgrClock = 0.0;
 }
 
 if ( mgrInited )
   {
    mgrTimer.tmAddr = NULL;
    mgrTimer.tmCount = 0;
    mgrTimer.tmWakeUp = 0;
    mgrTimer.tmReserved = 0;
    InsTime( (QElemPtr)&mgrTimer );
    PrimeTime( (QElemPtr)&mgrTimer, -MAX_TIME );
   }
 return( mgrClock );
}
#endif

/***********************************************************/
/*  Parsytec GCel timer.                                   */
/*  Provided by: Georg Wambach, gw@informatik.uni-koeln.de */
/***********************************************************/
#ifdef PARIX
#include <sys/time.h>

double dtime()
{
 double q;

 q = (double) (TimeNowHigh()) / (double) CLK_TCK_HIGH;

 return q;
}
#endif

/************************************************/
/*  Sun Solaris POSIX dtime() routine           */
/*  Provided by: Case Larsen, CTLarsen.lbl.gov  */
/************************************************/
#ifdef POSIX
#include <sys/time.h>
#include <sys/resource.h>
#include <sys/rusage.h>

#ifdef __hpux
#include <sys/syscall.h>
#endif

struct rusage rusage;

double dtime()
{
 double q;

 getrusage(RUSAGE_SELF,&rusage);

 q = (double)(rusage.ru_utime.tv_sec);
 q = q + (double)(rusage.ru_utime.tv_nsec) * 1.0e-09;
	
 return q;
}
#endif


/****************************************************/
/*  Windows NT (32 bit) dtime() routine             */
/*  Provided by: Piers Haken, piersh@microsoft.com  */
/****************************************************/
#ifdef WIN32
#include <windows.h>

double dtime(void)
{
 double q;

 q = (double)GetTickCount() * 1.0e-03;
	
 return q;
}
#endif

/*****************************************************/
/* Time according to POSIX.1  -  <J.Pelan@qub.ac.uk> */
/* Ref: "POSIX Programmer's Guide"  O'Reilly & Assoc.*/
/*****************************************************/
#ifdef POSIX1
#define _POSIX_SOURCE 1
#include <unistd.h>
#include <limits.h>
#include <sys/times.h>

struct tms tms;

double dtime()
{
 double q;
 times(&tms);
 q = (double)tms.tms_utime / (double)CLK_TCK;
 return q;
}
#endif

/* ---------- End of shuffle.c, Say goodnight Carrie ---------- */