// 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 8static 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 - 1unsigned 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] );       }}