Main Page | Namespace List | Class Hierarchy | Class List | File List | Class Members | File Members

Sobol.cpp

Go to the documentation of this file.
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 }

Note: Generated nightly - reload for latest version
Generated on Thu Dec 22 23:12:37 2005 for terreneuve by doxygen 1.3.6