quasic.sh01

Nistnet网络仿真软件

---- Cut Here and unpack ----
#!/bin/sh
# shar:	Shell Archiver  (v1.22)
#
#                                                                          
#                                                                          
#
#	Run the following text with /bin/sh to create:
#	  quasi.h
#	  quasi.c
#
if test -r s2_seq_.tmp
then echo "Must unpack archives in sequence!"
     next=`cat s2_seq_.tmp`; echo "Please unpack part $next next"
     exit 1; fi
echo "x - extracting quasi.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > quasi.h &&
X/* quasi.h		Quasi-Random Number generator, an object oriented
X                        implementation in C.
X
X		returns an n-dimensional vector of values in 0.0..1.0
X		maximum n is currently hardwired to 52
X
X			See W.H. Press and S.A. Teukolsky, 1989, Quasi- (that is,
X			Sub-) Random Numbers, Computers in Physics V3, No. 6,
X			(Nov/Dec 1989), pp. 76-79
X
X
X                rcsid: @(#)quasi.h	1.5 10:15:46 4/18/94   EFC
X
X*/
X
X#ifndef QUASI_RANDOM_H_
X#define QUASI_RANDOM_H_ 1.5
X
X
Xtypedef struct       /* data structure to maintain internal information */
X{                    /* this allows multiple independent generators to exist in
X                        a single application (each one would have its own
X			QRStruct) */
X
X      int err_flag;             /* err_flag == 0 if all is well */
X      int dim;
X      unsigned long int index;
X      unsigned long int *ix;
X} QRStruct;
X
X#ifdef __cplusplus
Xextern "C" {
X#endif
X
X/* initialized the data structure for this dimension */
Xint QuasiRandomInitialize(QRStruct* qr, int dimension);
X
X/* release internal storage for this QRStruct */
Xvoid QuasiRandomRelease(QRStruct* qr);
X
X/* get an n-dimensional quasi-random number */
Xvoid QuasiRandomNumber(QRStruct* qr, float* x);
X
X
X#ifdef __cplusplus
X}
X#endif
X
X#endif
X
X
X
X
X
SHAR_EOF
chmod 0644 quasi.h || echo "restore of quasi.h fails"
echo "x - extracting quasi.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > quasi.c &&
X/* quasi.c		Implementation of the Quasi-Random Number generator
X			currently hardwired to no more than 52 dimensions
X
X			See W.H. Press and S.A. Teukolsky, 1989, Quasi- (that is,
X			Sub-) Random Numbers, Computers in Physics V3, No. 6,
X			(Nov/Dec 1989), pp. 76-79
X
X
X*/
X
Xstatic char rcsid[] = "@(#)quasi.c	1.5 10:15:59 4/18/94   EFC";
X
X#include <stdio.h>
X
X#include "quasi.h"
X
X/* the primitive polynomial coefficients for up to degree 8  */
Xstatic int ip[] = { 0, 1, 1, 2, 1, 4, 2, 4, 7, 11, 13, 14,
X			  1, 13, 16, 19, 22, 25,
X			  1, 4, 7, 8, 14, 19, 21, 28, 31, 32, 37, 41, 42,
X			  50, 55, 56, 59, 62,
X			  14, 21, 22, 38, 47, 49, 50, 52, 56, 67, 70, 84,
X			  97, 103, 115, 122 };
Xstatic int mdeg[] = { 1, 2, 3, 3, 4, 4, 5, 5, 5, 5, 5, 5,
X			    6, 6, 6, 6, 6, 6,
X			    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
X			    8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8
X			     };
X
Xstatic int            maxdim = sizeof(mdeg) / sizeof(int);
X
Xstatic int            maxbit = 30;	    /* must be no more than
X					       number of bits in ulong - 1 */
X
Xstatic unsigned long  int*  iv = NULL;
Xstatic double               factor = 1.0;
Xstatic int instances = 0;
X
X#define INDEX(k,j)	[(k) + (j-1) * maxdim]
X
X
Xstatic void QuasiInit()	               /* initialize the direction numbers */
X{                                      /* only done once */
X	int j, k, l, ipp, niv;
X	unsigned long int mval, limit;
X	unsigned long int i;
X        
X        iv = (unsigned long int *)malloc( (niv = maxdim * maxbit)
X					 * sizeof( unsigned long int) );
X
X	for (k = 0; k < niv; k++)
X        	iv[k] = 0;
X
X        for (k = 0; k < maxdim; k++)
X        	iv[k] = 1;
X
X	mval = 4;
X	ipp = 1;
X
X        limit = (1L << (maxbit-1) );
X
X        for (k = maxdim, j = 0; k < niv-1; k += 2)
X        {
X		iv[k] = ipp;
X                if (++j == maxdim)
X                {
X		        if ( mval < limit )
X			          mval *= 2;
X                        ipp += 2;
X                        j = 0;
X                }
X
X		if ( ipp > mval )
X                	ipp = 1;
X                        
X                iv[k+1] = ipp;
X                if (++j == maxdim)
X                {
X		        if ( mval < limit )
X			          mval *= 2;
X                        ipp += 2;
X                        j = 0;
X                }
X		else
X                {
X                	ipp += 2;
X                        if ( ipp > mval )
X                        	ipp = 1;
X                }
X                
X
X        }
X
X	for (k = 0; k < maxdim; k++)
X        {
X        	/* normalize the set iv values */
X        	for (j = 1; j <= mdeg[k]; j++)
X                	iv INDEX(k,j) *= (1L << (maxbit - j));
X
X		/* calcululate the rest of the iv values */
X		for (j = mdeg[k] + 1; j <= maxbit; j++)
X                {
X                	ipp = ip[k];
X                        
X                        /* calculate Gray code of iv */
X                        i = iv INDEX(k, j - mdeg[k]);
X                        i ^= i / (1L << mdeg[k]);
X                        
X                        for (l = mdeg[k] - 1; l >= 1; l--)
X                        {
X                        	if ( ipp & 1 )
X                                	i ^= iv INDEX(k, j-l);
X                                ipp /= 2;
X                        }
X
X                        iv INDEX(k,j) = i;
X
X        	}
X	}
X
X        factor = 1.0 / (1L << maxbit);
X
X                        
X}
X
Xint QuasiRandomInitialize(QRStruct* qr, int dimension)
X{
X        int k;
X
X        qr->dim = dimension;
X	qr->index = 0;
X
X	if ( qr->dim > maxdim )		/* if dimension is too large   */
X        {				/* truncate and set error flag */
X        	qr->dim = maxdim;
X		qr->err_flag = -1;
X        }
X        else
X        	qr->err_flag = 0;
X                                
X	qr->ix = (unsigned long int *)malloc(qr->dim * sizeof( unsigned long int) );
X
X	for (k = 0; k < qr->dim; k++)
X        	qr->ix[k] = 0L;
X
X                
X	if ( instances++ == 0 )
X		QuasiInit();
X
X	return qr->err_flag;
X}
X
Xvoid QuasiRandomRelease(QRStruct* qr)
X{
X	if ( --instances == 0 )
X	    free( iv );
X
X	free( qr->ix );
X
X	qr->err_flag = -9;
X}
X
Xvoid QuasiRandomNumber(QRStruct* qr, float* x)
X{
X	int i, j, k;
X	unsigned long int im = qr->index++;
X
X        /* find rightmost zero bit  */
X        for (j = 0; j < maxbit; j++, im >>= 1)
X        	if ( (im & 1L) == 0 )
X                		break;
X
X	i = j * maxdim;
X
X        for (k = 0; k < qr->dim; k++)
X        {
X        	qr->ix[k] ^= iv[i + k];		/* integer values  */
X                x[k]   = (float) ( factor * (double)qr->ix[k] );
X        }
X
X}
X
X
X#ifdef MAIN
X
X/* Demonstration test code, writes quasirandom numbers to standard out
X
X
X   To compile with the test driver:
X
X   cc -DMAIN -o quasi quasi.c -lm
X
X
X
X   usage:
X
X      quasi [points [dimension]]
X
X
X*/
X
Xint main(int argc, char** argv)
X{
X    int i, k;
X    float *x;
X    QRStruct qr;
X    int dim = 0;
X    int points = 0;
X
X    if ( argc > 1 )
X          points = atoi( argv[1] );
X
X    if ( argc > 2 )
X          dim = atoi( argv[2] );
X
X    if ( points < 1 )
X           points = 100;    /* default to 100 points */
X
X    if ( dim < 1 )
X           dim = 2;         /* default to 2 dimensional points */
X
X
X    x = (float *)malloc( dim * sizeof( float ) );
X
X    QuasiRandomInitialize( &qr, dim );
X
X    for (i = 0; i < points; i++)
X    {
X              QuasiRandomNumber( &qr, x );
X              for (k = 0; k < dim; k++)
X                       printf("%f  ", x[k]);
X	      putchar( '\n' );
X    }
X
X
X    QuasiRandomRelease( &qr );
X
X    free( x );
X
X}
X
X#endif
X
X
X
SHAR_EOF
chmod 0644 quasi.c || echo "restore of quasi.c fails"
exit 0