Randomness of MSC 5.1 rand()
Pat T Shea
pats at cup.portal.com
Tue May 9 20:28:47 AEST 1989
'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;
}
More information about the Comp.lang.c
mailing list