HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_SobolSequence.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_SobolSequence.h ( UT Library, C++)
7  *
8  * COMMENTS: A sobol low-discrepancy sequence.
9  */
10 
11 #ifndef __UT_SobolSequence__
12 #define __UT_SobolSequence__
13 
14 #include "UT_API.h"
15 #include "UT_Array.h"
16 #include "UT_Assert.h"
17 #include "UT_StackAlloc.h"
18 #include "UT_Vector3.h"
19 #include <SYS/SYS_Math.h>
20 #include <SYS/SYS_StaticAssert.h>
21 #include <SYS/SYS_Types.h>
22 #include <type_traits>
23 
24 /// Random Sobol sequence of dimension _DIM_ with _BITS_ bits of integral
25 /// *unsigned* type _itype_
26 template <class itype, uint DIM, uint BITS = sizeof(itype)*8>
28 {
29 public:
31 
32  constexpr uint dimension() const { return DIM; }
33  constexpr uint bits() const { return BITS; }
34 
35  /// Advance to the next sample, storing the result into x. When myN
36  /// exceeds the range of the sequence, it will automatically wrap.
37  void advance();
38 
39  /// Return the sequence value for dimension dim.
40  itype x(exint dim) const { return myX[dim]; }
41  fpreal32 xf(exint dim) const
42  {
43  SYS_FPRealUnionF fval;
44 
45  // Move significant digit to line up with mantissa and
46  // bitwise OR it with 1.0 to get [1.0,2.0).
47  // BITS_FROM_23 is moved out to avoid compiler warning.
48  if (BITS < 23)
49  fval.uval = 0x3f800000 | (uint)(x(dim) << BITS_FROM_23);
50  else
51  fval.uval = 0x3f800000 | (uint)(x(dim) >> BITS_FROM_23);
52  return fval.fval - 1;
53  }
54 
55  /// Reset the sequence
56  void rewind(itype n = 0);
57 
58 private:
59  /// Returns an odd random value between 0 and val.
60  static constexpr uint
61  orand(uint val, uint &seed)
62  {
63  return ((SYSfastRandomInt(seed) % (val >> 1)) << 1) | 1;
64  }
65  static constexpr itype
66  computeMask()
67  {
68  itype mask = 0;
69  for (uint i = 0; i < BITS; i++)
70  mask |= itype(1) << i;
71  return mask;
72  }
73 
74 private:
75  static constexpr uint LEN_IV = DIM * (BITS + 1);
76  static constexpr uint LEN_X = DIM;
77  static constexpr uint MAXDIM = 6;
78  static constexpr itype MASK = computeMask();
79  static constexpr uint BITS_FROM_23 = (BITS<23) ? (23-BITS) : (BITS-23);
80 
81  UT_Array<itype> myIV;
82  UT_Array<itype> myX;
83  itype myN;
84 };
85 
87 {
89 
90 public:
91  UT_VectorSequence() = default;
92 
93  using PARENT::advance;
94  using PARENT::rewind;
95 
97  {
98  return UT_Vector3(PARENT::xf(0),
99  PARENT::xf(1),
100  PARENT::xf(2));
101  }
102 };
103 
104 
105 //////////////////////////////////////////////////////////////////////////////
106 //
107 // Implementation
108 //
109 
110 template <class itype, uint DIM, uint BITS>
112  : myIV(LEN_IV, LEN_IV)
113  , myX(LEN_X, LEN_X)
114  , myN(0)
115 {
116  constexpr uint mdeg[MAXDIM] = {1,2,3,3,4,4};
117  constexpr uint ip[MAXDIM] = {0,1,1,2,1,4};
118 
120  SYS_STATIC_ASSERT(DIM >= 1 && DIM <= MAXDIM);
121  SYS_STATIC_ASSERT(BITS >= 1 && BITS <= sizeof(itype)*8);
122 
123  // Indexing
124  itype **iu = (itype **)UTstackAlloc((BITS+1)*sizeof(itype *));
125  for (uint i = 0; i < BITS+1; i++)
126  iu[i] = &myIV[i*DIM];
127 
128  // Fill in the initial odd random values
129  uint seed = 1337;
130  for (uint j = 0; j < DIM; j++)
131  for (uint i = 0; i < mdeg[j]; i++)
132  iu[i][j] = orand(1 << (i+1), seed);
133 
134  for (uint k = 0; k < DIM; k++)
135  {
136  // Initial values only require normalization
137  for (uint j = 0; j < mdeg[k]; j++)
138  iu[j][k] <<= (BITS-j-1);
139 
140  // Calculate the remaining values using the recurrence equation
141  for (uint j = mdeg[k]; j < BITS+1; j++)
142  {
143  itype ipp = ip[k];
144 
145  // First and last terms
146  itype tmp = iu[j-mdeg[k]][k];
147  tmp ^= (tmp >> mdeg[k]);
148 
149  // Intermediate terms
150  for (uint l = mdeg[k]-1; l >= 1; l--)
151  {
152  if (ipp & 1)
153  tmp ^= iu[j-l][k];
154  ipp >>= 1;
155  }
156  iu[j][k] = tmp;
157  }
158  }
159 
160  rewind();
161 
162  UTstackFree(iu);
163 }
164 
165 template <class itype, uint DIM, uint BITS>
166 void
168 {
169  if (myN < MASK)
170  {
171  // Find the first 0-bit in n
172  uint i;
173  itype n = myN & MASK;
174  for (i = 0; i < BITS; i++)
175  {
176  if (!(n & 1))
177  break;
178  n >>= 1;
179  }
180 
181  // XOR the bits into the sequence entry
182  n = i*DIM;
183  for (i = 0; i < DIM; i++)
184  {
185  myX[i] ^= myIV[n+i];
186  }
187  myN++;
188  }
189  else
190  {
191  rewind();
192  }
193 }
194 
195 template <class itype, uint DIM, uint BITS>
196 void
198 {
199  n &= MASK;
200  myN = n;
201  for (uint i = 0; i < DIM; i++)
202  myX[i] = 0;
203 
204  uint off = 0;
205  while (n)
206  {
207  if ((n+1) & 0x2)
208  {
209  for (uint i = 0; i < DIM; i++)
210  myX[i] ^= myIV[off+i];
211  }
212  n >>= 1;
213  off += DIM;
214  }
215 }
216 
217 #endif // __UT_SobolSequence__
#define SYS_STATIC_ASSERT(expr)
GLsizei const GLfloat * value
Definition: glcorearb.h:824
UT_Vector3T< float > UT_Vector3
int64 exint
Definition: SYS_Types.h:125
float fpreal32
Definition: SYS_Types.h:200
UT_VectorSequence()=default
GLdouble GLdouble x2
Definition: glad.h:2349
#define UTstackAlloc(size)
Definition: UT_StackAlloc.h:34
GLdouble n
Definition: glcorearb.h:2008
fpreal32 xf(exint dim) const
constexpr uint bits() const
GLint GLuint mask
Definition: glcorearb.h:124
UT_Vector3 getVector()
GLint j
Definition: glad.h:2733
void rewind(itype n=0)
Reset the sequence.
GLuint GLfloat * val
Definition: glcorearb.h:1608
itype x(exint dim) const
Return the sequence value for dimension dim.
unsigned int uint
Definition: SYS_Types.h:45
constexpr uint dimension() const
#define UTstackFree(ptr)
Definition: UT_StackAlloc.h:38