2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2017,2018, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * \brief Implementation of the 2x64 ThreeFry random engine
39 * \author Erik Lindahl <erik.lindahl@gmail.com>
41 * \ingroup module_random
44 #ifndef GMX_RANDOM_THREEFRY_H
45 #define GMX_RANDOM_THREEFRY_H
50 #include "gromacs/math/functions.h"
51 #include "gromacs/random/seed.h"
52 #include "gromacs/utility/classhelpers.h"
53 #include "gromacs/utility/exceptions.h"
56 * The GROMACS implementation of the ThreeFry random engine has been
57 * heavily inspired by the versions proposed to Boost by:
59 * John Salmon, Copyright 2010-2014 by D. E. Shaw Research
60 * https://github.com/DEShawResearch/Random123-Boost
62 * Thijs van den Berg, Copyright (c) 2014 M.A. (Thijs) van den Berg
63 * https://github.com/sitmo/threefry
65 * Both of them are covered by the Boost Software License:
67 * Boost Software License - Version 1.0 - August 17th, 2003
69 * Permission is hereby granted, free of charge, to any person or organization
70 * obtaining a copy of the software and accompanying documentation covered by
71 * this license (the "Software") to use, reproduce, display, distribute,
72 * execute, and transmit the Software, and to prepare derivative works of the
73 * Software, and to permit third-parties to whom the Software is furnished to
74 * do so, all subject to the following:
76 * The copyright notices in the Software and this entire statement, including
77 * the above license grant, this restriction and the following disclaimer,
78 * must be included in all copies of the Software, in whole or in part, and
79 * all derivative works of the Software, unless such copies or derivative
80 * works are solely in the form of machine-executable object code generated by
81 * a source language processor.
83 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
84 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
85 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
86 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
87 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
88 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
89 * DEALINGS IN THE SOFTWARE.
97 // Variable-bitfield counters used to increment internal counters as
98 // part of std::arrays.
103 /*! \brief Clear highBits higest bits of ctr, return false if they were non-zero.
105 * This function clears the space required for the internal counters,
106 * and returns true if they were correctly zero when calling, false otherwise.
108 * \tparam UIntType Integer type to use for each word in counter
109 * \tparam words Number of UIntType words in counter
110 * \tparam highBits Number of bits to check. The template parameter makes it
111 * possible to optimize this extensively at compile time.
112 * \param ctr Reference to counter to check and clear.
114 template<class UIntType, std::size_t words, unsigned int highBits>
116 checkAndClear(std::array<UIntType, words> * ctr)
118 const std::size_t bitsPerWord = std::numeric_limits<UIntType>::digits;
119 const std::size_t bitsTotal = bitsPerWord*words;
121 static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");
123 const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord;
124 const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
125 const UIntType lastWordOne = static_cast<UIntType>(1) << lastWordLowBitIdx;
126 const UIntType mask = lastWordOne-1;
130 for (unsigned int i = words-1; i > lastWordIdx; --i)
138 if (highBits > 0 && (*ctr)[lastWordIdx] >= lastWordOne)
141 (*ctr)[lastWordIdx] &= mask;
146 /*! \brief Increment the internal counter in highBits by one
148 * \tparam UIntType Integer type to use for each word in counter
149 * \tparam words Number of UIntType words in counter
150 * \tparam highBits Number of bits reserved for the internal counter.
151 * \param ctr Reference to the counter value to increment.
153 * \throws InternalError if internal counter space is exhausted.
155 * This routine will work across the word boundaries for any number
156 * of internal counter bits that fits in the total counter.
158 template<class UIntType, std::size_t words, unsigned int highBits>
160 increment(std::array<UIntType, words> * ctr)
162 const std::size_t bitsPerWord = std::numeric_limits<UIntType>::digits;
163 const std::size_t bitsTotal = bitsPerWord*words;
165 static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");
167 const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord;
168 const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
169 const UIntType lastWordOne = static_cast<UIntType>(1) << lastWordLowBitIdx;
171 // For algorithm & efficiency reasons we need to store the internal counter in
172 // the same array as the user-provided counter, so we use the higest bits, possibly
173 // crossing several words.
175 // To have the computer help us with the dirty carry arithmetics we store the bits
176 // in the internal counter part in normal fashion, but the internal counter words in
177 // reverse order; the highest word of the total counter array (words-1) is thus
178 // the least significant part of the internal counter (if it spans several words).
180 // The incrementation works as follows:
182 // 0) If the index of the least significant internal counter word is larger
183 // than words-1, there was never any space.
184 // 1) If the internal counter spans more than one word, we must have one or
185 // more internal counter words that correspond entirely to the this counter.
186 // Start with the least significant one (words-1) and increment it.
187 // If the new value is not zero we did not loop around (no carry), so everything
188 // is good, and we are done - return!
189 // If the new value is zero, we need to move the carry result to the next word,
190 // so we just continue the loop until we have gone through all words that
191 // are internal-counter-only.
192 // 2) After the loop, there is stuff remaining to add, and by definition there
193 // is some internal counter space in the next word, but the question
194 // is if we have exhausted it. We already created a constant that corresponds
195 // to the bit that represents '1' for the internal counter part of this word.
196 // When we add this constant it will not affect the user-counter-part at all,
197 // and if we exhaust the internal counter space the high bits will cause the entire
198 // word to wrap around, and the result will be smaller than the bit we added.
199 // If this happens we throw, otherwise we're done.
201 // Since all constants will be evaluated at compile time, this entire routine
202 // will usually be reduced to simply incrementing a word by a constant, and throwing
203 // if the result is smaller than the constant.
205 if (lastWordIdx >= words)
207 GMX_THROW(InternalError("Cannot increment random engine defined with 0 internal counter bits."));
210 for (unsigned int i = words-1; i > lastWordIdx; --i)
215 return; // No carry means we are done
218 (*ctr)[lastWordIdx] += lastWordOne;
219 if ((*ctr)[lastWordIdx] < lastWordOne)
221 GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
225 /*! \brief Increment the internal counter in highBits by a value.
227 * \tparam UIntType Integer type to use for each word in counter
228 * \tparam words Number of UIntType words in counter
229 * \tparam highBits Number of bits reserved for the internal counter.
230 * \param ctr Reference to the counter to increment.
231 * \param addend Value to add to internal.
233 * \throws InternalError if internal counter space is exhausted.
235 * This routine will work across the word boundaries for any number
236 * of internal counter bits that fits in the total counter.
238 template<class UIntType, std::size_t words, unsigned int highBits>
240 increment(std::array<UIntType, words> * ctr, UIntType addend)
242 const std::size_t bitsPerWord = std::numeric_limits<UIntType>::digits;
243 const std::size_t bitsTotal = bitsPerWord*words;
245 static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");
247 const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord;
248 const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
249 const UIntType lastWordOne = static_cast<UIntType>(1) << lastWordLowBitIdx;
250 const UIntType lastWordMaxVal = (~static_cast<UIntType>(0)) >> lastWordLowBitIdx;
252 if (lastWordIdx >= words)
254 GMX_THROW(InternalError("Cannot increment random engine defined with 0 internal counter bits."));
257 for (unsigned int i = words-1; i > lastWordIdx; --i)
260 addend = ((*ctr)[i] < addend); // 1 is the carry!
267 if (addend > lastWordMaxVal)
269 GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
271 addend *= lastWordOne;
273 (*ctr)[lastWordIdx] += addend;
275 if ((*ctr)[lastWordIdx] < addend)
277 GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
283 /*! \brief General implementation class for ThreeFry counter-based random engines.
285 * This class is used to implement several different ThreeFry2x64 random engines
286 * differing in the number of rounds executed in and the number of bits reserved
287 * for the internal counter. It is compatible with C++11 random engines, and
288 * can be used e.g. with all random distributions from the standard library.
290 * ThreeFry is a counter-based rather than state-based random engine. This
291 * means that we seed it with a "key", after which we can get the
292 * N:th random number in a sequence (specified by a counter) directly. This
293 * means we are guaranteed the same sequence of numbers even when running in
294 * parallel if using e.g. step and atom index as counters.
296 * However, it is also useful to be able to use it as a normal random engine,
297 * for instance if you need more than 2 64-bit random values for a specific
298 * counter value, not to mention where you just need good normal random numbers.
299 * To achieve this, this implementation uses John Salmon's idea of reserving
300 * a couple of the highest bits in the user-provided counter for an internal
301 * counter. For instance, if reserving 3 bits, this means you get a stream of
302 * 8 iterations (each with 2 random values) after every restart. If you call
303 * the engine after these bits have been exhausted, it will throw an
304 * exception to make sure you don't get overlapping streams by mistake.
305 * Reserving 3 bits also means you can only use 64-3=61 bits of the highest
306 * word when restarting (i.e., setting) the counters.
308 * This version also supports using internalCounterBits=0. In this case the
309 * random engine will be able to return a single counter round, i.e. 2 64-bit
310 * values for ThreeFry2x64, after which an exception is thrown. In this case no
311 * high bits are reserved, which means the class implements the raw ThreeFry2x64
314 * \tparam rounds The number of encryption iterations used when generating.
315 * This can in principle be any value, but 20 rounds has been
316 * shown to pass all BigCrush random tests, and with 13 rounds
317 * only one fails. This is a very stringent test, and the
318 * standard Mersenne Twister engine fails two, so 13 rounds
319 * should be a perfectly fine balance in most cases.
320 * \tparam internalCounterBits
321 * Number of high bits in the user-provided counter reserved
322 * for the internal counter. The number of values the engine
323 * can return after each restart will be
324 * words*2^internalCounterBits.
326 template<unsigned int rounds, unsigned int internalCounterBits>
327 class ThreeFry2x64General
329 // While this class will formally work with any value for rounds, there is
330 // no reason to go lower than 13, and this might help catch some typos.
331 // If we find a reason to use lower values in the future, or if you simply
332 // want to test, this assert can safely be removed.
333 static_assert(rounds >= 13, "You should not use less than 13 encryption rounds for ThreeFry2x64.");
336 // result_type must be lower case to be compatible with C++11 standard library
338 /*! \brief Integer type for output. */
339 typedef uint64_t result_type;
340 /*! \brief Use array for counter & key states so it is allocated on the stack */
341 typedef std::array<result_type, 2> counter_type;
345 /*! \brief Rotate value left by specified number of bits
347 * \param i Value to rotate (result_type, which should be 64-bit).
348 * \param bits Number of bits to rotate i.
350 * \return Input value rotated 'bits' left.
353 rotLeft(result_type i, unsigned int bits)
355 return (i << bits) | (i >> (std::numeric_limits<result_type>::digits-bits));
358 /*! \brief Perform encryption step for ThreeFry2x64 algorithm
360 * It performs the encryption step of the standard ThreeFish symmetric-key
361 * tweakable block cipher, which is the core of the ThreeFry random
362 * engine. The number of encryption rounds is specified by the class
363 * template parameter 'rounds'.
365 * \param key Reference to key value
366 * \param ctr Counter value to use
368 * \return Newly encrypted 2x64 block, according to the class template parameters.
371 generateBlock(const counter_type &key,
372 const counter_type &ctr)
374 const unsigned int rotations[] = {16, 42, 12, 31, 16, 32, 24, 21};
375 counter_type x = ctr;
377 result_type ks[3] = { 0x0, 0x0, 0x1bd11bdaa9fc1a22 };
379 // This is actually a pretty simple routine that merely executes the
380 // for-block specified further down 'rounds' times. However, both
381 // clang and gcc have problems unrolling and replacing rotations[r%8]
382 // with constants, so we unroll the first 20 iterations manually.
386 ks[0] = key[0]; ks[2] ^= key[0]; x[0] = x[0] + key[0];
387 ks[1] = key[1]; ks[2] ^= key[1]; x[1] = x[1] + key[1];
388 x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0];
390 if (rounds > 1) { x[0] += x[1]; x[1] = rotLeft(x[1], 42); x[1] ^= x[0]; }
391 if (rounds > 2) { x[0] += x[1]; x[1] = rotLeft(x[1], 12); x[1] ^= x[0]; }
392 if (rounds > 3) { x[0] += x[1]; x[1] = rotLeft(x[1], 31); x[1] ^= x[0]; x[0] += ks[1]; x[1] += ks[2] + 1; }
393 if (rounds > 4) { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
394 if (rounds > 5) { x[0] += x[1]; x[1] = rotLeft(x[1], 32); x[1] ^= x[0]; }
395 if (rounds > 6) { x[0] += x[1]; x[1] = rotLeft(x[1], 24); x[1] ^= x[0]; }
396 if (rounds > 7) { x[0] += x[1]; x[1] = rotLeft(x[1], 21); x[1] ^= x[0]; x[0] += ks[2]; x[1] += ks[0] + 2; }
397 if (rounds > 8) { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
398 if (rounds > 9) { x[0] += x[1]; x[1] = rotLeft(x[1], 42); x[1] ^= x[0]; }
399 if (rounds > 10) { x[0] += x[1]; x[1] = rotLeft(x[1], 12); x[1] ^= x[0]; }
400 if (rounds > 11) { x[0] += x[1]; x[1] = rotLeft(x[1], 31); x[1] ^= x[0]; x[0] += ks[0]; x[1] += ks[1] + 3; }
401 if (rounds > 12) { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
402 if (rounds > 13) { x[0] += x[1]; x[1] = rotLeft(x[1], 32); x[1] ^= x[0]; }
403 if (rounds > 14) { x[0] += x[1]; x[1] = rotLeft(x[1], 24); x[1] ^= x[0]; }
404 if (rounds > 15) { x[0] += x[1]; x[1] = rotLeft(x[1], 21); x[1] ^= x[0]; x[0] += ks[1]; x[1] += ks[2] + 4; }
405 if (rounds > 16) { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
406 if (rounds > 17) { x[0] += x[1]; x[1] = rotLeft(x[1], 42); x[1] ^= x[0]; }
407 if (rounds > 18) { x[0] += x[1]; x[1] = rotLeft(x[1], 12); x[1] ^= x[0]; }
408 if (rounds > 19) { x[0] += x[1]; x[1] = rotLeft(x[1], 31); x[1] ^= x[0]; x[0] += ks[2]; x[1] += ks[0] + 5; }
410 for (unsigned int r = 20; r < rounds; r++)
413 x[1] = rotLeft(x[1], rotations[r%8]);
415 if (( (r + 1) & 3 ) == 0)
417 unsigned int r4 = (r + 1) >> 2;
418 x[0] += ks[ r4 % 3 ];
419 x[1] += ks[ (r4 + 1) % 3 ] + r4;
426 //! \brief Smallest value that can be returned from random engine.
427 #if !defined(_MSC_VER)
430 // Avoid constexpr bug in MSVC 2015, note that max() below does work
433 result_type min() { return std::numeric_limits<result_type>::min(); }
435 //! \brief Largest value that can be returned from random engine.
437 result_type max() { return std::numeric_limits<result_type>::max(); }
439 /*! \brief Construct random engine with 2x64 key values
441 * This constructor takes two values, and should only be used with
442 * the 2x64 implementations.
444 * \param key0 Random seed in the form of a 64-bit unsigned value.
445 * \param domain Random domain. This is used to guarantee that different
446 * applications of a random engine inside the code get different
447 * streams of random numbers, without requiring the user
448 * to provide lots of random seeds. Pick a value from the
449 * RandomDomain class, or RandomDomain::Other if it is
450 * not important. In the latter case you might want to use
451 * \ref gmx::DefaultRandomEngine instead.
453 * \note The random domain is really another 64-bit seed value.
455 * \throws InternalError if the high bits needed to encode the number of counter
458 ThreeFry2x64General(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other)
463 /*! \brief Construct random engine from 2x64-bit unsigned integers
465 * This constructor assigns the raw 128 bit key data from unsigned integers.
466 * It is meant for the case when you want full control over the key,
467 * for instance to compare with reference values of the ThreeFry
468 * function during testing.
470 * \param key0 First word of key/random seed.
471 * \param key1 Second word of key/random seed.
473 * \throws InternalError if the high bits needed to encode the number of counter
474 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
476 ThreeFry2x64General(uint64_t key0, uint64_t key1)
481 /*! \brief Seed 2x64 random engine with two 64-bit key values
483 * \param key0 First word of random seed, in the form of 64-bit unsigned values.
484 * \param domain Random domain. This is used to guarantee that different
485 * applications of a random engine inside the code get different
486 * streams of random numbers, without requiring the user
487 * to provide lots of random seeds. Pick a value from the
488 * RandomDomain class, or RandomDomain::Other if it is
489 * not important. In the latter case you might want to use
490 * \ref gmx::DefaultRandomEngine instead.
492 * \note The random domain is really another 64-bit seed value.
494 * Re-initialized the seed similar to the counter constructor.
495 * Same rules apply: The highest few bits of the last word are
496 * reserved to encode the number of internal counter bits, but
497 * to save the user the trouble of making sure these are zero
498 * when using e.g. a random device, we just ignore them.
501 seed(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other)
503 seed(key0, static_cast<uint64_t>(domain));
506 /*! \brief Seed random engine from 2x64-bit unsigned integers
508 * This assigns the raw 128 bit key data from unsigned integers.
509 * It is meant for the case when you want full control over the key,
510 * for instance to compare with reference values of the ThreeFry
511 * function during testing.
513 * \param key0 First word of key/random seed.
514 * \param key1 Second word of key/random seed.
516 * \throws InternalError if the high bits needed to encode the number of counter
517 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
520 seed(uint64_t key0, uint64_t key1)
522 const unsigned int internalCounterBitsBits = (internalCounterBits > 0) ? ( StaticLog2<internalCounterBits>::value + 1 ) : 0;
524 key_ = {{key0, key1}};
526 if (internalCounterBits > 0)
528 internal::highBitCounter::checkAndClear<result_type, 2, internalCounterBitsBits>(&key_);
529 internal::highBitCounter::increment<result_type, 2, internalCounterBitsBits>(&key_, internalCounterBits-1);
534 /*! \brief Restart 2x64 random engine counter from 2 64-bit values
536 * \param ctr0 First word of new counter, in the form of 64-bit unsigned values.
537 * \param ctr1 Second word of new counter
539 * Restarting the engine with a new counter is extremely fast with ThreeFry64,
540 * and basically just consists of storing the counter value, so you should
541 * use this liberally in your innermost loops to restart the engine with
542 * e.g. the current step and atom index as counter values.
544 * \throws InternalError if any of the highest bits that are reserved
545 * for the internal part of the counter are set. The number of
546 * reserved bits is to the last template parameter to the class.
549 restart(uint64_t ctr0 = 0, uint64_t ctr1 = 0)
552 counter_ = {{ctr0, ctr1}};
553 if (!internal::highBitCounter::checkAndClear<result_type, 2, internalCounterBits>(&counter_))
555 GMX_THROW(InternalError("High bits of counter are reserved for the internal stream counter."));
557 block_ = generateBlock(key_, counter_);
561 /*! \brief Generate the next random number
563 * This will return the next stored 64-bit value if one is available,
564 * and otherwise generate a new block, update the internal counters, and
565 * return the first value while storing the others.
567 * \throws InternalError if the internal counter space is exhausted.
572 if (index_ >= c_resultsPerCounter_)
574 internal::highBitCounter::increment<result_type, 2, internalCounterBits>(&counter_);
575 block_ = generateBlock(key_, counter_);
578 return block_[index_++];
581 /*! \brief Skip next n random numbers
583 * Moves the internal random stream for the give key/counter value
584 * n positions forward. The count is based on the number of random values
585 * returned, such that skipping 5 values gives exactly the same result as
586 * drawing 5 values that are ignored.
588 * \param n Number of values to jump forward.
590 * \throws InternalError if the internal counter space is exhausted.
595 index_ += n % c_resultsPerCounter_;
596 n /= c_resultsPerCounter_;
598 if (index_ > c_resultsPerCounter_)
600 index_ -= c_resultsPerCounter_;
604 // Make sure the state is the same as if we came to this counter and
605 // index by natural generation.
606 if (index_ == 0 && n > 0)
608 index_ = c_resultsPerCounter_;
611 internal::highBitCounter::increment<result_type, 2, internalCounterBits>(&counter_, n);
612 block_ = generateBlock(key_, counter_);
615 /*! \brief Return true if two ThreeFry2x64 engines are identical
617 * \param x Instance to compare with.
619 * This routine should return true if the two engines will generate
620 * identical random streams when drawing.
623 operator==(const ThreeFry2x64General<rounds, internalCounterBits> &x) const
625 // block_ is uniquely specified by key_ and counter_.
626 return (key_ == x.key_ && counter_ == x.counter_ && index_ == x.index_);
629 /*! \brief Return true of two ThreeFry2x64 engines are not identical
631 * \param x Instance to compare with.
633 * This routine should return true if the two engines will generate
634 * different random streams when drawing.
637 operator!=(const ThreeFry2x64General<rounds, internalCounterBits> &x) const { return !operator==(x); }
641 /*! \brief Number of results returned for each invocation of the block generation */
642 static const unsigned int c_resultsPerCounter_ = static_cast<unsigned int>(sizeof(counter_type)/sizeof(result_type));
644 /*! \brief ThreeFry2x64 key, i.e. the random seed for this stream.
646 * The highest few bits of the key are replaced to encode the value of
647 * internalCounterBits, in order to make all streams unique.
651 /*! \brief ThreeFry2x64 total counter.
653 * The highest internalCounterBits are reserved for an internal counter
654 * so that the combination of a key and counter provides a stream that
655 * returns 2*2^internalCounterBits (ThreeFry2x64) random 64-bit values before
656 * the internal counter space is exhausted and an exception is thrown.
658 counter_type counter_;
659 /*! \brief The present block encrypted from values of key and counter. */
661 /*! \brief Index of the next value in block_ to return from random engine */
664 GMX_DISALLOW_COPY_AND_ASSIGN(ThreeFry2x64General);
668 /*! \brief ThreeFry2x64 random engine with 20 iteractions.
670 * \tparam internalCounterBits, default 64.
672 * This class provides very high quality random numbers that pass all
673 * BigCrush tests, it works with two 64-bit values each for keys and
674 * counters, and is most efficient when we only need a few random values
675 * before restarting the counters with new values.
677 template<unsigned int internalCounterBits = 64>
678 class ThreeFry2x64 : public ThreeFry2x64General<20, internalCounterBits>
681 /*! \brief Construct ThreeFry random engine with 2x64 key values, 20 rounds.
683 * \param key0 Random seed in the form of a 64-bit unsigned value.
684 * \param domain Random domain. This is used to guarantee that different
685 * applications of a random engine inside the code get different
686 * streams of random numbers, without requiring the user
687 * to provide lots of random seeds. Pick a value from the
688 * RandomDomain class, or RandomDomain::Other if it is
689 * not important. In the latter case you might want to use
690 * \ref gmx::DefaultRandomEngine instead.
692 * \note The random domain is really another 64-bit seed value.
694 * \throws InternalError if the high bits needed to encode the number of counter
697 ThreeFry2x64(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) : ThreeFry2x64General<20, internalCounterBits>(key0, domain) {}
699 /*! \brief Construct random engine from 2x64-bit unsigned integers, 20 rounds
701 * This constructor assigns the raw 128 bit key data from unsigned integers.
702 * It is meant for the case when you want full control over the key,
703 * for instance to compare with reference values of the ThreeFry
704 * function during testing.
706 * \param key0 First word of key/random seed.
707 * \param key1 Second word of key/random seed.
709 * \throws InternalError if the high bits needed to encode the number of counter
710 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
712 ThreeFry2x64(uint64_t key0, uint64_t key1) : ThreeFry2x64General<20, internalCounterBits>(key0, key1) {}
715 /*! \brief ThreeFry2x64 random engine with 13 iteractions.
717 * \tparam internalCounterBits, default 64.
719 * This class provides relatively high quality random numbers that only
720 * fail one BigCrush test, and it is a bit faster than the 20-round version.
721 * It works with two 64-bit values each for keys and counters, and is most
722 * efficient when we only need a few random values before restarting
723 * the counters with new values.
725 template<unsigned int internalCounterBits = 64>
726 class ThreeFry2x64Fast : public ThreeFry2x64General<13, internalCounterBits>
729 /*! \brief Construct ThreeFry random engine with 2x64 key values, 13 rounds.
731 * \param key0 Random seed in the form of a 64-bit unsigned value.
732 * \param domain Random domain. This is used to guarantee that different
733 * applications of a random engine inside the code get different
734 * streams of random numbers, without requiring the user
735 * to provide lots of random seeds. Pick a value from the
736 * RandomDomain class, or RandomDomain::Other if it is
737 * not important. In the latter case you might want to use
738 * \ref gmx::DefaultRandomEngine instead.
740 * \note The random domain is really another 64-bit seed value.
742 * \throws InternalError if the high bits needed to encode the number of counter
745 ThreeFry2x64Fast(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) : ThreeFry2x64General<13, internalCounterBits>(key0, domain) {}
747 /*! \brief Construct ThreeFry random engine from 2x64-bit unsigned integers, 13 rounds.
749 * This constructor assigns the raw 128 bit key data from unsigned integers.
750 * It is meant for the case when you want full control over the key,
751 * for instance to compare with reference values of the ThreeFry
752 * function during testing.
754 * \param key0 First word of key/random seed.
755 * \param key1 Second word of key/random seed.
757 * \throws InternalError if the high bits needed to encode the number of counter
758 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
760 ThreeFry2x64Fast(uint64_t key0, uint64_t key1) : ThreeFry2x64General<13, internalCounterBits>(key0, key1) {}
765 /*! \brief Default fast and accurate random engine in Gromacs
767 * This engine will return 2*2^64 random results using the default
768 * gmx::RandomDomain::Other stream, and can be initialized with a single
769 * seed argument without having to remember empty template angle brackets.
771 typedef ThreeFry2x64Fast<> DefaultRandomEngine;
775 #endif // GMX_RANDOM_THREEFRY_H