00001 #include ".\sobol.h" 00002 #include <minmax.h> 00003 00004 Sobol::Sobol(LongNatural Seed_) 00005 { 00006 Seed=Seed_; 00007 n_=-1; 00008 sobseq(&n_,x_); 00009 n_=1; 00010 } 00011 00012 Sobol::~Sobol() 00013 { 00014 } 00015 00016 void Sobol::SetSeed(LongNatural Seed_) { 00017 Seed=Seed_; 00018 n_=-1; 00019 sobseq(&n_,x_); 00020 n_=1; 00021 } 00022 00023 void Sobol::sobseq(Integer *n, Real x[]) 00024 { 00025 Natural j,k,l; 00026 LongInteger i,im,ipp; 00027 static LongNatural mdeg[MAXDIM+1]={0,1,2,3,3,4,4}; 00028 static LongNatural in; 00029 static LongNatural ix[MAXDIM+1]; 00030 static LongNatural *iu[MAXBIT+1]; 00031 static LongNatural ip[MAXDIM+1]={0,0,1,1,2,1,4}; 00032 static LongNatural iv[MAXDIM*MAXBIT+1]= 00033 {0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9}; 00034 static Real fac; 00035 00036 if (*n < 0) { 00037 for (k=1;k<MAXDIM;k++) ix[k]=0; 00038 in=0; 00039 if (iv[1] != 1) return; 00040 fac=1.0/(1L << MAXBIT); 00041 for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k]; 00042 for (k=1;k<=MAXDIM;k++) { 00043 for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j); 00044 for (j=mdeg[k]+1;j<=MAXBIT;j++) { 00045 ipp=ip[k]; 00046 i=iu[j-mdeg[k]][k]; 00047 i ^= (i >> mdeg[k]); 00048 for (l=mdeg[k]-1;l>=1;l--) { 00049 if (ipp & 1) i ^= iu[j-l][k]; 00050 ipp >>= 1; 00051 } 00052 iu[j][k]=i; 00053 } 00054 } 00055 } else { 00056 im=in++; 00057 for (j=1;j<=MAXBIT;j++) { 00058 if (!(im & 1)) break; 00059 im >>= 1; 00060 } 00061 if (j > MAXBIT) 00062 j=MAXBIT; 00063 im=(j-1)*MAXDIM; 00064 for (k=1;k<=min(*n,MAXDIM);k++) { 00065 ix[k] ^= iv[im+k]; 00066 x[k]=ix[k]*fac; 00067 } 00068 } 00069 } 00070 00071 VeryLongNatural Sobol::Max() { 00072 return 2; 00073 } 00074 00075 LongNatural Sobol::Min() { 00076 return 1; 00077 } 00078 00079 LongNatural Sobol::GetOneRandomInteger() 00080 { 00081 return 1; 00082 } 00083 00084 Real Sobol::getUniform() 00085 { 00086 sobseq(&n_,x_); 00087 return x_[1]; 00088 }