Fermat
lfsr.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2010-2016, NVIDIA Corporation
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  * * Redistributions of source code must retain the above copyright
8  * notice, this list of conditions and the following disclaimer.
9  * * Redistributions in binary form must reproduce the above copyright
10  * notice, this list of conditions and the following disclaimer in the
11  * documentation and/or other materials provided with the distribution.
12  * * Neither the name of NVIDIA Corporation nor the
13  * names of its contributors may be used to endorse or promote products
14  * derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  */
27 
34 #pragma once
35 
36 #include <cugar/basic/types.h>
37 
38 
39 namespace cugar {
40 
43 
55 
66 class LFSRGeneratorMatrix // Markov chain quasi-Monte Carlo linear feedback shift register.
67 {
68 public:
69  enum Offset_type
70  {
71  ORIGINAL, // The original ones from the paper.
72  GOOD_PROJECTIONS // Maximized minimum distance for some projections.
73  };
74 
75  // Construct a small LFSR with period 2^m - 1, for 3 <= m <= 32.
76  // The offset_type describes which set of offset values to use.
77  CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
79 
80  // Construct a small LFSR with period 2^m - 1, for 3 <= m <= 32.
81  // The offset_type describes which set of offset values to use.
82  CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
83  explicit LFSRGeneratorMatrix(uint32 m, Offset_type offset_type = GOOD_PROJECTIONS);
84 
85  // Update the given state, and return the next value in the sequence,
86  // using the scramble value as a random bit shift. For the first invocation
87  // an arbitrary state value 0 < *state < 2^m may be used. The scramble parameter
88  // may be an arbitrary (possibly random) value. To generate the scrambled
89  // coordinates of the origin it is valid to pass *state == 0, but in this case
90  // the state will not be updated.
91  CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
92  float next(uint32 scramble, uint32* state) const;
93 
94 public:
95  uint32 m_m; // period 2^m - 1
96  uint32 m_f[32]; // f_d^s
97 };
98 
110 {
111  // Set the initial state, using the scramble value as a random bit shift. For the first invocation
112  // an arbitrary state value 0 < *state < 2^m may be used. The scramble parameter
113  // may be an arbitrary (possibly random) value. To generate the scrambled
114  // coordinates of the origin it is valid to pass *state == 0, but in this case
115  // the state will not be updated.
116  CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
117  LFSRRandomStream(LFSRGeneratorMatrix* matrix, const uint32 state, const uint32 scramble = 1361u) :
118  m_matrix(matrix), m_state(state ? state : 0xFFFFFFFF), m_scramble(scramble) {}
119 
120  CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
121  float next();
122 
123  LFSRGeneratorMatrix* m_matrix;
124  uint32 m_state;
125  uint32 m_scramble;
126 };
127 
128 CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
129 LFSRGeneratorMatrix::LFSRGeneratorMatrix(const uint32 m, const LFSRGeneratorMatrix::Offset_type offset_type)
130  : m_m(m)
131 {
132  // Table taken from T. Hansen and G. Mullen:
133  // "Primitive Polynomials over Finite Fields".
134  // It is implied that the coefficient for t^m is 1.
135  const uint32 primitive_polynomials[32 - 3 + 1] =
136  {
137  (1 << 1) | 1, // 3
138  (1 << 1) | 1, // 4
139  (1 << 2) | 1, // 5
140  (1 << 1) | 1, // 6
141  (1 << 1) | 1, // 7
142  (1 << 4) | (1 << 3) | (1 << 2) | 1, // 8
143  (1 << 4) | 1, // 9
144  (1 << 3) | 1, // 10
145  (1 << 2) | 1, // 11
146  (1 << 6) | (1 << 4) | (1 << 1) | 1, // 12
147  (1 << 4) | (1 << 3) | (1 << 1) | 1, // 13
148  (1 << 5) | (1 << 3) | (1 << 1) | 1, // 14
149  (1 << 1) | 1, // 15
150  (1 << 5) | (1 << 3) | (1 << 2) | 1, // 16
151  (1 << 3) | 1, // 17
152  (1 << 7) | 1, // 18
153  (1 << 5) | (1 << 2) | (1 << 1) | 1, // 19
154  (1 << 3) | 1, // 20
155  (1 << 2) | 1, // 21
156  (1 << 1) | 1, // 22
157  (1 << 5) | 1, // 23
158  (1 << 4) | (1 << 3) | (1 << 1) | 1, // 24
159  (1 << 3) | 1, // 25
160  (1 << 6) | (1 << 2) | (1 << 1) | 1, // 26
161  (1 << 5) | (1 << 2) | (1 << 1) | 1, // 27
162  (1 << 3) | 1, // 28
163  (1 << 2) | 1, // 29
164  (1 << 6) | (1 << 4) | (1 << 1) | 1, // 30
165  (1 << 3) | 1, // 31
166  (1 << 7) | (1 << 6) | (1 << 2) | 1 // 32
167  };
168 
169  // The original offsets 10 <= m <= 32 are taken from S. Chen, M. Matsumoto, T. Nishimura,
170  // and A. Owen: "New inputs and methods for Markov chain quasi-Monte Carlo".
171  // The alternative set of offsets for 3 <= m <= 16 was computed by Leonhard Gruenschloss.
172  // They should also yield maximal equidistribution as described in P. L'Ecuyer: "Maximally
173  // Equidistributed Combined Tausworthe Generators", but their projections might have better
174  // maximized minimum distance properties.
175  const uint32 offsets[32 - 3 + 1][2] =
176  {
177  // org / good proj.
178  { 1, 1 }, // 3
179  { 2, 2 }, // 4
180  { 15, 15 }, // 5
181  { 8, 8 }, // 6
182  { 4, 4 }, // 7
183  { 41, 41 }, // 8
184  { 113, 113 }, // 9
185  { 115, 226 }, // 10 *
186  { 291, 520 }, // 11 *
187  { 172, 1583 }, // 12 *
188  { 267, 2242 }, // 13 *
189  { 332, 2312 }, // 14 *
190  { 388, 38 }, // 15 *
191  { 283, 13981 }, // 16 *
192  { 514, 514 }, // 17
193  { 698, 698 }, // 18
194  { 706, 706 }, // 19
195  { 1304, 1304 }, // 20
196  { 920, 920 }, // 21
197  { 1336, 1336 }, // 22
198  { 1236, 1236 }, // 23
199  { 1511, 1511 }, // 24
200  { 1445, 1445 }, // 25
201  { 1906, 1906 }, // 26
202  { 1875, 1875 }, // 27
203  { 2573, 2573 }, // 28
204  { 2633, 2633 }, // 29
205  { 2423, 2423 }, // 30
206  { 3573, 3573 }, // 31
207  { 3632, 3632 } // 32
208  };
209 
210  // Construct the matrix that corresponds to a single transition.
211  uint32 matrix[32];
212  matrix[m - 1] = 0;
213  for (uint32 i = 1, pp = primitive_polynomials[m - 3]; i < m; ++i, pp >>= 1)
214  {
215  matrix[m - 1] |= (pp & 1) << (m - i); // Reverse bits.
216  matrix[i - 1] = 1 << (m - i - 1);
217  }
218 
219  // Apply the matrix exponentiation according to the offset.
220  uint32 result0[32], result1[32]; // Storage for temporary results.
221  for (unsigned i = 0; i < m; ++i)
222  result0[i] = matrix[i]; // Copy over row.
223  uint32* in = result0;
224  uint32* out = result1;
225  const uint32 offset = offsets[m - 3][static_cast<int>(offset_type)];
226  for (unsigned i = 1; i < offset; ++i)
227  {
228  // Perform matrix multiplication: out = in * matrix.
229  for (uint32 y = 0; y < m; ++y)
230  {
231  out[y] = 0;
232  for (uint32 x = 0; x < m; ++x)
233  for (uint32 i = 0; i < m; ++i)
234  out[y] ^= (((in[y] >> i) & (matrix[m - i - 1] >> x)) & 1) << x;
235  }
236 
237  // Swap input and output.
238  unsigned* tmp = in;
239  in = out;
240  out = tmp;
241  }
242 
243  // Transpose the result for simpler multiplication.
244  for (uint32 y = 0; y < m; ++y)
245  {
246  m_f[y] = 0;
247  for (uint32 x = 0; x < m; ++x)
248  m_f[y] |= ((in[x] >> y) & 1) << (m - x - 1);
249  }
250 }
251 
252 CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
253 float LFSRGeneratorMatrix::next(const uint32 scramble, uint32* const state) const
254 {
255  //check(state);
256 
257  // Compute the matrix-vector multiplication using one matrix column at a time.
258  uint32 result = 0;
259  for (uint32 i = 0, s = *state; s; ++i, s >>= 1)
260  {
261  if (s & 1)
262  result ^= m_f[i];
263  }
264 
265 #if !defined(FLT_EPSILON)
266  const float FLT_EPSILON = 1.0e-10f;
267 #endif
268 
269  //check(result <= ~(1u << m_m));
270  *state = result; // Write result back.
271  result = (result << (32 - m_m)) ^ scramble; // Apply scrambling.
272  const float fresult = result * (1.f / float(uint64(1ULL) << 32)); // Map to [0, 1).
273  return fresult <= 1.0f - FLT_EPSILON ?
274  fresult : 1.0f - FLT_EPSILON;
275 }
276 
277 CUGAR_HOST_DEVICE CUGAR_FORCEINLINE
278 float LFSRRandomStream::next()
279 {
280  return m_matrix->next(m_scramble, &m_state);
281 }
282 
285 
286 } // namespace cugar
Definition: lfsr.h:66
Definition: lfsr.h:109
Define a vector_view POD type and plain_view() for std::vector.
Definition: diff.h:38