Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!ames!hc!lll-winken!uunet!portal!cup.portal.com!pats
From: pats@cup.portal.com (Pat T Shea)
Newsgroups: comp.lang.c
Subject: Re: Randomness of MSC 5.1 rand()
Message-ID: <18148@cup.portal.com>
Date: 9 May 89 10:28:47 GMT
References: <1598@oregon.uoregon.edu>
Organization: The Portal System (TM)
Lines: 203

'asked the same question and hung the following Quick 'n Dirty together
. . . . .   MSC v5.1
 
/* ----------SnipSnipSnipSnipSnipSnip--------- */
 
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys\types.h>
#include <sys\timeb.h>
 
#define MAX_RAND 100
 
extern void main( int argc, char * *argv );
extern int random( int r) ;
extern void randinit( void );
extern void moment( double data[], register int n, double *ave,
                    double *adev, double *sdev, double *svar,
                    double *skew, double *curt );
 
void main( int argc, char **argv )
{
   register int i, j;
   const int k = atoi( argv[1] );
   unsigned long r[MAX_RAND];
   double hits[MAX_RAND + 1], average = 0.0, *p_average, ave_dev = 0.0,
          *p_ave_dev, std_dev = 0.0, *p_std_dev,
          std_var = 0.0, *p_std_var, skew = 0.0, *p_skew,
          kurtosis = 0.0, *p_kurtosis;
   time_t start, end;
 
   if ( argc < 2 || k < 2 || k > MAX_RAND ) {
      printf( "\a\n\tSYNTAX:  t_rand <n>\n\n\t" );
      printf( "Any number, greater than 1 but less " );
      printf( "than or equal to %d !!.....\n\n", MAX_RAND );
      exit( 1 );
   }
 
/* Initialize the Random Number generator based on SysTime            */
 
   randinit();
 
/* Get start time                                                     */
 
   time( &start );
 
/* Set pointers and zero out accumulator buckets....                  */
 
   p_average = &average;
   p_ave_dev = &ave_dev;
   p_std_dev = &std_dev;
   p_std_var = &std_var;
   p_skew = &skew;
   p_kurtosis = &kurtosis;
   for ( j = 0; j < MAX_RAND; j++ ) {
      r[j] = 0;
   }
   for ( j = 0; j < MAX_RAND + 1; j++ ) {
      hits[j] = 0.0;
   }
 
/* and loop like a mad....                                            */
 
   printf( "\n\t\tRunning 500,000 iterations....\n\n" );
   printf( "\n\tThis takes 'bout 5 minutes on an IBM XT @ 4.77MHz\n" );
   printf( "\t\tand 17-18 seconds on a COMPAQ 386/20 !!\n\n" );
   for ( j = 0; j < 20; j++ ) {
      for ( i = 0; i < 25000; i++ ) {
         r[random( k )]++;
      }
   }
   printf( "\aFor %d Random Numbers in the range, 0 to %d .....\n\n\t",
            k, k - 1 );
   for ( j = 0; j > MAX_RAND; j++ ) {
      if ( r[j] ) {
         printf( "%4d was hit...... %9lu times or %7.4f%%.\n\t",
                  j, r[j], ( (double)( r[j] ) / 5000 ));
      }
   }
 
/* Load up 'double' array and calculate statistics                    */
 
   for ( j = 0; j < MAX_RAND; j++ ) {
      hits[j + 1] = (double) r[j];
   }
   moment( hits, k, p_average, p_ave_dev, p_std_dev, p_std_var,
           p_skew, p_kurtosis );
   printf( "\tMEAN of Hits was %13.5f\n", average );
   printf( "\t\tAVE DEV of Hits was %10.2f\n", ave_dev );
   printf( "\t\tSTD DEV of Hits was %10.2f\n", std_dev );
   printf( "\t\tVARIANCE of Hits was %9.2f\n", std_var );
   printf( "\tSKEWNESS was %f, ", skew );
   printf( "and KURTOSIS was %f.\n", kurtosis );
   printf( "\tThe STD DEV as a Percent of the MEAN is %7.5f%%.\n\n",
            100 * std_dev / average );
 
   time( &end );
   printf( "\t\tThe Elapsed Time was %.0f seconds.\n\n", difftime( end, start ));
   exit( 0 );
}
 
/***********************************************************************
 *
 *                   to generate random integers
 *                in the range of 0 to ( numb - 1 ).......
 *                    /Pat Shea @ Psi! 11-19-87
 *
 *****/
 
int random( register int numb )
{
   long int raw_rand;
 
   raw_rand = (long) rand();
   return( (int)( raw_rand * numb / 32768 ));
}
 
/***********************************************************************
 *
 *                         void randinit( void )
 *
 *    Initializes the Random Number generator with an unsigned integer
 *       in the range 0 to 65535.  This 'seed' number is derived by
 *       checking the System Clock, determining the number of SECONDS
 *       since January 01, 1970, dividing this seconds number by 65535,
 *       and 'seeding' with the remainder from the division.  OR.....
 *       there are 65,000+ possible 'seeds' or initializers for the
 *       random sequences thus generated.  Other techniques involve
 *       seeding off the thousandths of a second on the System Clock
 *       but this is NOT bullet proof in that not all systems maintain
 *       that parameter, AND there are only 1000 possible seeds.  If
 *       that paramenter is not present, the random number generator
 *       will always present the same sequence, seeded of the base
 *       of '1'.  This risk in using that parameter is TOO HIGH.
 *
 *                    /Pat Shea @ Psi! 11-19-87
 *****/
 
void randinit( void )
{
   struct timeb sys_time;
 
   ftime( &sys_time );
   srand( (unsigned int) ( sys_time.time % 65535L ));
   return;
}
 
/***********************************************************************
 *
 *  NAME:
 *    moment.c - crunches elementary statistics
 *
 *    adapted
 *       from  - Numeric Recipes in C - The Art of Scientific Computing
 *               by Press, Flannery, Teukolsky and Vetterling
 *               Cambridge University Press 1988  p475f.
 *
 *             - bulletpruf'd and ANSI'd /PatS 06-15-88
 *
 *  SYNOPSIS:
 *       void moment( double data[], register int n, double *ave,
 *                    double *adev, double *sdev, double *svar,
 *                    double *skew, double *curt )
 *
 ****/
 
#include <conio.h>
#include <math.h>
 
void moment( double data[], register int n, double *ave, double *adev,
             double *sdev, double *svar, double *skew, double *curt )
{
   register int j;
   double p, s = 0.0;
 
   *adev = ( *svar ) = ( *skew ) = ( *curt ) = 0.0;
   if ( n <= 1 ) {
      (void) cputs( "\n\n\rMUST have at least 2 " );
      (void) cputs( "observations for \"MOMENT\" to work !!\n\n\r" );
   }
   for ( j = 1; j <= n; j++ ) {
      s += data[j];
   }
   *ave = ( s / n );
   for ( j = 1; j <= n; j++ ) {
      *adev += fabs( s = data[j] - ( *ave ));
      *svar += ( p = s * s );
      *skew += ( p *= s );
      *curt += ( p *= s );
   }
   *adev /= n;
   *svar /= ( n - 1 );
   *sdev = sqrt( *svar );
   if ( *svar ) {
      *skew /= ( n * ( *svar ) * ( *sdev ));
      *curt = ( *curt ) / ( n * ( *svar ) * ( *svar )) - 3.0;
   }
   else {
      (void) cputs( "\n\n\rCannot calculate skew/kurtosis when " );
      (void) cputs( "variance is 0.0 !! ( in \"MOMENT\" )\n\n\r" );
   }
   return;
}