quasi.cpp

伪随机数生成器

// quasi.cpp		Implementation of the Quasi-Random Number generator
//			currently hardwired to no more than 52 dimensions

//                      (c) Copyright 1996, Everett F. Carter Jr.
//                      Permission is granted by the author to use
//                      this software for any application provided this
//                      copyright notice is preserved.

static const char rcsid[] = "@(#)quasi.c++	1.6 10:55:50 4/3/96   EFC";

#ifdef DEBUG
#include <iostream.h>
#endif
#include <stdlib.h>

#include <quasi.hpp>

// the primitive polynomial coefficients for up to degree 8
static const int ip[] = { 0, 1, 1, 2, 1, 4, 2, 4, 7, 11, 13, 14,
			  1, 13, 16, 19, 22, 25,
			  1, 4, 7, 8, 14, 19, 21, 28, 31, 32, 37, 41, 42,
			  50, 55, 56, 59, 62,
			  14, 21, 22, 38, 47, 49, 50, 52, 56, 67, 70, 84,
			  97, 103, 115, 122 };
static const int mdeg[] = { 1, 2, 3, 3, 4, 4, 5, 5, 5, 5, 5, 5,
			    6, 6, 6, 6, 6, 6,
			    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
			    8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8
			     };

const int            QuasiRandom::maxdim = sizeof(mdeg) / sizeof(int);

const int            QuasiRandom::maxbit = 30;	    // must be no more than
						    // number of bits in ulong - 1
unsigned long  int*  QuasiRandom::iv = NULL;
double               QuasiRandom::factor = 1.0;
int                  QuasiRandom::instances = 0;

#define INDEX(k,j)	[(k) + (j-1) * maxdim]


void QuasiRandom::init()	// initialize the direction numbers
{
	int j, k, l, ipp, niv;
	unsigned long int i, mval;
        
        iv = new unsigned long int[niv = maxdim * maxbit];

	for (k = 0; k < niv; k++)
        	iv[k] = 0;

        for (k = 0; k < maxdim; k++)
        	iv[k] = 1;

	mval = 4;
	ipp = 1;
        
        for (k = maxdim, j = 0; k < niv-1; k += 2)
        {
		iv[k] = ipp;
                if (++j == maxdim)
                {
			mval *= 2;
                        ipp += 2;
                        j = 0;
                }

		if ( ipp > mval )
                	ipp = 1;
                        
                iv[k+1] = ipp;
                if (++j == maxdim)
                {
			mval *= 2;
                        ipp += 2;
                        j = 0;
                }
		else
                {
                	ipp += 2;
                        if ( ipp > mval )
                        	ipp = 1;
                }
                

        }

	for (k = 0; k < maxdim; k++)
        {
        	// normalize the set iv values
        	for (j = 1; j <= mdeg[k]; j++)
                	iv INDEX(k,j) *= (1L << (maxbit - j));

		// calcululate the rest of the iv values
		for (j = mdeg[k] + 1; j <= maxbit; j++)
                {
                	ipp = ip[k];
                        
                        // calculate Gray code of iv
                        i = iv INDEX(k, j - mdeg[k]);
                        i ^= i / (1L << mdeg[k]);
                        
                        for (l = mdeg[k] - 1; l >= 1; l--)
                        {
                        	if ( ipp & 1 )
                                	i ^= iv INDEX(k, j-l);
                                ipp /= 2;
                        }

                        iv INDEX(k,j) = i;

        	}
	}

        factor = 1.0 / (1L << maxbit);

                        
}

QuasiRandom::QuasiRandom(const int dimension) : dim(dimension), index(0)
{

	if ( dim > maxdim )		// if dimension is too large
        {				// truncate and set error flag
        	dim = maxdim;
		err_flag = -1;
        }
        else
        	err_flag = 0;
                                
	ix = new unsigned long int[dim];

	for (int k = 0; k < dim; k++)
        	ix[k] = 0L;

                
	if ( instances++ == 0 )
		init();
                
}

QuasiRandom::~QuasiRandom()
{
	if ( --instances == 0 )
#if defined( __ZTC__ ) && __ZTC__ <= 0x301
        	delete [maxbit*maxdim]iv;
#else
		delete []iv;
#endif

#if defined( __ZTC__ ) && __ZTC__ <= 0x301
        delete [dim]ix;
#else
	delete []ix;
#endif
        
}

void QuasiRandom::number(float* x)
{
	int i, j, k;
	unsigned long int im = index++;

        // find rightmost zero bit
        for (j = 0; j < maxbit; j++, im >>= 1)
        	if ( (im & 1L) == 0 )
                		break;

	i = j * maxdim;

        for (k = 0; k < dim; k++)
        {
        	ix[k] ^= iv[i + k];		// integer values
                x[k]   = (float) ( factor * (double)ix[k] );
       }

}

void QuasiRandom::number(BasicArray& x)
{
	int i, j, k;
	unsigned long int im = index++;

        // find rightmost zero bit
        for (j = 0; j < maxbit; j++, im >>= 1)
        	if ( (im & 1L) == 0 )
                		break;

	i = j * maxdim;

	x.resize( dim );

        for (k = 0; k < dim; k++)
        {
        	ix[k] ^= iv[i + k];		// integer values
                x[k]   = (float) ( factor * (double)ix[k] );
       }

}