2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2017,2018,2019,2021, 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
51 #include "gromacs/math/functions.h"
52 #include "gromacs/random/seed.h"
53 #include "gromacs/utility/classhelpers.h"
54 #include "gromacs/utility/exceptions.h"
57 * The GROMACS implementation of the ThreeFry random engine has been
58 * heavily inspired by the versions proposed to Boost by:
60 * John Salmon, Copyright 2010-2014 by D. E. Shaw Research
61 * https://github.com/DEShawResearch/Random123-Boost
63 * Thijs van den Berg, Copyright (c) 2014 M.A. (Thijs) van den Berg
64 * https://github.com/sitmo/threefry
66 * Both of them are covered by the Boost Software License:
68 * Boost Software License - Version 1.0 - August 17th, 2003
70 * Permission is hereby granted, free of charge, to any person or organization
71 * obtaining a copy of the software and accompanying documentation covered by
72 * this license (the "Software") to use, reproduce, display, distribute,
73 * execute, and transmit the Software, and to prepare derivative works of the
74 * Software, and to permit third-parties to whom the Software is furnished to
75 * do so, all subject to the following:
77 * The copyright notices in the Software and this entire statement, including
78 * the above license grant, this restriction and the following disclaimer,
79 * must be included in all copies of the Software, in whole or in part, and
80 * all derivative works of the Software, unless such copies or derivative
81 * works are solely in the form of machine-executable object code generated by
82 * a source language processor.
84 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
85 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
86 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
87 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
88 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
89 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
90 * DEALINGS IN THE SOFTWARE.
98 // Variable-bitfield counters used to increment internal counters as
99 // part of std::arrays.
101 struct highBitCounter
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>
115 static bool checkAndClear(std::array<UIntType, words>* ctr)
117 const std::size_t bitsPerWord = std::numeric_limits<UIntType>::digits;
118 const std::size_t bitsTotal = bitsPerWord * words;
120 static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");
122 const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord;
123 const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
124 const UIntType lastWordOne = static_cast<UIntType>(1) << lastWordLowBitIdx;
125 const UIntType mask = lastWordOne - 1;
129 for (unsigned int i = words - 1; i > lastWordIdx; --i)
137 if (highBits > 0 && (*ctr)[lastWordIdx] >= lastWordOne)
140 (*ctr)[lastWordIdx] &= mask;
145 /*! \brief Increment the internal counter in highBits by one
147 * \tparam UIntType Integer type to use for each word in counter
148 * \tparam words Number of UIntType words in counter
149 * \tparam highBits Number of bits reserved for the internal counter.
150 * \param ctr Reference to the counter value to increment.
152 * \throws InternalError if internal counter space is exhausted.
154 * This routine will work across the word boundaries for any number
155 * of internal counter bits that fits in the total counter.
157 template<class UIntType, std::size_t words, unsigned int highBits>
158 static void increment(std::array<UIntType, words>* ctr)
160 const std::size_t bitsPerWord = std::numeric_limits<UIntType>::digits;
161 const std::size_t bitsTotal = bitsPerWord * words;
163 static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");
165 const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord;
166 const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
167 const UIntType lastWordOne = static_cast<UIntType>(1) << lastWordLowBitIdx;
169 // For algorithm & efficiency reasons we need to store the internal counter in
170 // the same array as the user-provided counter, so we use the higest bits, possibly
171 // crossing several words.
173 // To have the computer help us with the dirty carry arithmetics we store the bits
174 // in the internal counter part in normal fashion, but the internal counter words in
175 // reverse order; the highest word of the total counter array (words-1) is thus
176 // the least significant part of the internal counter (if it spans several words).
178 // The incrementation works as follows:
180 // 0) If the index of the least significant internal counter word is larger
181 // than words-1, there was never any space.
182 // 1) If the internal counter spans more than one word, we must have one or
183 // more internal counter words that correspond entirely to the this counter.
184 // Start with the least significant one (words-1) and increment it.
185 // If the new value is not zero we did not loop around (no carry), so everything
186 // is good, and we are done - return!
187 // If the new value is zero, we need to move the carry result to the next word,
188 // so we just continue the loop until we have gone through all words that
189 // are internal-counter-only.
190 // 2) After the loop, there is stuff remaining to add, and by definition there
191 // is some internal counter space in the next word, but the question
192 // is if we have exhausted it. We already created a constant that corresponds
193 // to the bit that represents '1' for the internal counter part of this word.
194 // When we add this constant it will not affect the user-counter-part at all,
195 // and if we exhaust the internal counter space the high bits will cause the entire
196 // word to wrap around, and the result will be smaller than the bit we added.
197 // If this happens we throw, otherwise we're done.
199 // Since all constants will be evaluated at compile time, this entire routine
200 // will usually be reduced to simply incrementing a word by a constant, and throwing
201 // if the result is smaller than the constant.
203 if (lastWordIdx >= words)
205 GMX_THROW(InternalError(
206 "Cannot increment random engine defined with 0 internal counter bits."));
209 for (unsigned int i = words - 1; i > lastWordIdx; --i)
214 return; // No carry means we are done
217 (*ctr)[lastWordIdx] += lastWordOne;
218 if ((*ctr)[lastWordIdx] < lastWordOne)
220 GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
224 /*! \brief Increment the internal counter in highBits by a value.
226 * \tparam UIntType Integer type to use for each word in counter
227 * \tparam words Number of UIntType words in counter
228 * \tparam highBits Number of bits reserved for the internal counter.
229 * \param ctr Reference to the counter to increment.
230 * \param addend Value to add to internal.
232 * \throws InternalError if internal counter space is exhausted.
234 * This routine will work across the word boundaries for any number
235 * of internal counter bits that fits in the total counter.
237 template<class UIntType, std::size_t words, unsigned int highBits>
238 static void increment(std::array<UIntType, words>* ctr, UIntType addend)
240 const std::size_t bitsPerWord = std::numeric_limits<UIntType>::digits;
241 const std::size_t bitsTotal = bitsPerWord * words;
243 static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");
245 const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord;
246 const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
247 const UIntType lastWordOne = static_cast<UIntType>(1) << lastWordLowBitIdx;
248 const UIntType lastWordMaxVal = (~static_cast<UIntType>(0)) >> lastWordLowBitIdx;
250 if (lastWordIdx >= words)
252 GMX_THROW(InternalError(
253 "Cannot increment random engine defined with 0 internal counter bits."));
256 for (unsigned int i = words - 1; i > lastWordIdx; --i)
259 addend = ((*ctr)[i] < addend); // 1 is the carry!
266 if (addend > lastWordMaxVal)
268 GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
270 addend *= lastWordOne;
272 (*ctr)[lastWordIdx] += addend;
274 if ((*ctr)[lastWordIdx] < addend)
276 GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
280 } // namespace internal
282 /*! \brief General implementation class for ThreeFry counter-based random engines.
284 * This class is used to implement several different ThreeFry2x64 random engines
285 * differing in the number of rounds executed in and the number of bits reserved
286 * for the internal counter. It is compatible with C++11 random engines, and
287 * can be used e.g. with all random distributions from the standard library.
289 * ThreeFry is a counter-based rather than state-based random engine. This
290 * means that we seed it with a "key", after which we can get the
291 * N:th random number in a sequence (specified by a counter) directly. This
292 * means we are guaranteed the same sequence of numbers even when running in
293 * parallel if using e.g. step and atom index as counters.
295 * However, it is also useful to be able to use it as a normal random engine,
296 * for instance if you need more than 2 64-bit random values for a specific
297 * counter value, not to mention where you just need good normal random numbers.
298 * To achieve this, this implementation uses John Salmon's idea of reserving
299 * a couple of the highest bits in the user-provided counter for an internal
300 * counter. For instance, if reserving 3 bits, this means you get a stream of
301 * 8 iterations (each with 2 random values) after every restart. If you call
302 * the engine after these bits have been exhausted, it will throw an
303 * exception to make sure you don't get overlapping streams by mistake.
304 * Reserving 3 bits also means you can only use 64-3=61 bits of the highest
305 * word when restarting (i.e., setting) the counters.
307 * This version also supports using internalCounterBits=0. In this case the
308 * random engine will be able to return a single counter round, i.e. 2 64-bit
309 * values for ThreeFry2x64, after which an exception is thrown. In this case no
310 * high bits are reserved, which means the class implements the raw ThreeFry2x64
313 * \tparam rounds The number of encryption iterations used when generating.
314 * This can in principle be any value, but 20 rounds has been
315 * shown to pass all BigCrush random tests, and with 13 rounds
316 * only one fails. This is a very stringent test, and the
317 * standard Mersenne Twister engine fails two, so 13 rounds
318 * should be a perfectly fine balance in most cases.
319 * \tparam internalCounterBits
320 * Number of high bits in the user-provided counter reserved
321 * for the internal counter. The number of values the engine
322 * can return after each restart will be
323 * words*2^internalCounterBits.
325 template<unsigned int rounds, unsigned int internalCounterBits>
326 class ThreeFry2x64General
328 // While this class will formally work with any value for rounds, there is
329 // no reason to go lower than 13, and this might help catch some typos.
330 // If we find a reason to use lower values in the future, or if you simply
331 // want to test, this assert can safely be removed.
332 static_assert(rounds >= 13,
333 "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;
344 /*! \brief Rotate value left by specified number of bits
346 * \param i Value to rotate (result_type, which should be 64-bit).
347 * \param bits Number of bits to rotate i.
349 * \return Input value rotated 'bits' left.
351 result_type rotLeft(result_type i, unsigned int bits)
353 return (i << bits) | (i >> (std::numeric_limits<result_type>::digits - bits));
356 /*! \brief Perform encryption step for ThreeFry2x64 algorithm
358 * It performs the encryption step of the standard ThreeFish symmetric-key
359 * tweakable block cipher, which is the core of the ThreeFry random
360 * engine. The number of encryption rounds is specified by the class
361 * template parameter 'rounds'.
363 * \param key Reference to key value
364 * \param ctr Counter value to use
366 * \return Newly encrypted 2x64 block, according to the class template parameters.
368 counter_type generateBlock(const counter_type& key, const counter_type& ctr)
370 const unsigned int rotations[] = { 16, 42, 12, 31, 16, 32, 24, 21 };
371 counter_type x = ctr;
373 result_type ks[3] = { 0x0, 0x0, 0x1bd11bdaa9fc1a22 };
375 // This is actually a pretty simple routine that merely executes the
376 // for-block specified further down 'rounds' times. However, both
377 // clang and gcc have problems unrolling and replacing rotations[r%8]
378 // with constants, so we unroll the first 20 iterations manually.
384 x[0] = x[0] + key[0];
387 x[1] = x[1] + key[1];
389 x[1] = rotLeft(x[1], 16);
395 x[1] = rotLeft(x[1], 42);
401 x[1] = rotLeft(x[1], 12);
407 x[1] = rotLeft(x[1], 31);
415 x[1] = rotLeft(x[1], 16);
421 x[1] = rotLeft(x[1], 32);
427 x[1] = rotLeft(x[1], 24);
433 x[1] = rotLeft(x[1], 21);
441 x[1] = rotLeft(x[1], 16);
447 x[1] = rotLeft(x[1], 42);
453 x[1] = rotLeft(x[1], 12);
459 x[1] = rotLeft(x[1], 31);
467 x[1] = rotLeft(x[1], 16);
473 x[1] = rotLeft(x[1], 32);
479 x[1] = rotLeft(x[1], 24);
485 x[1] = rotLeft(x[1], 21);
493 x[1] = rotLeft(x[1], 16);
499 x[1] = rotLeft(x[1], 42);
505 x[1] = rotLeft(x[1], 12);
511 x[1] = rotLeft(x[1], 31);
517 for (unsigned int r = 20; r < rounds; r++)
520 x[1] = rotLeft(x[1], rotations[r % 8]);
522 if (((r + 1) & 3) == 0)
524 unsigned int r4 = (r + 1) >> 2;
526 x[1] += ks[(r4 + 1) % 3] + r4;
533 //! \brief Smallest value that can be returned from random engine.
534 #if !defined(_MSC_VER)
537 // Avoid constexpr bug in MSVC 2015, note that max() below does work
543 return std::numeric_limits<result_type>::min();
546 //! \brief Largest value that can be returned from random engine.
547 static constexpr result_type max() { return std::numeric_limits<result_type>::max(); }
549 /*! \brief Construct random engine with 2x64 key values
551 * This constructor takes two values, and should only be used with
552 * the 2x64 implementations.
554 * \param key0 Random seed in the form of a 64-bit unsigned value.
555 * \param domain Random domain. This is used to guarantee that different
556 * applications of a random engine inside the code get different
557 * streams of random numbers, without requiring the user
558 * to provide lots of random seeds. Pick a value from the
559 * RandomDomain class, or RandomDomain::Other if it is
560 * not important. In the latter case you might want to use
561 * \ref gmx::DefaultRandomEngine instead.
563 * \note The random domain is really another 64-bit seed value.
565 * \throws InternalError if the high bits needed to encode the number of counter
568 //NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
569 ThreeFry2x64General(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other)
574 /*! \brief Construct random engine from 2x64-bit unsigned integers
576 * This constructor assigns the raw 128 bit key data from unsigned integers.
577 * It is meant for the case when you want full control over the key,
578 * for instance to compare with reference values of the ThreeFry
579 * function during testing.
581 * \param key0 First word of key/random seed.
582 * \param key1 Second word of key/random seed.
584 * \throws InternalError if the high bits needed to encode the number of counter
585 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
587 //NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init)
588 ThreeFry2x64General(uint64_t key0, uint64_t key1) { seed(key0, key1); }
590 /*! \brief Seed 2x64 random engine with two 64-bit key values
592 * \param key0 First word of random seed, in the form of 64-bit unsigned values.
593 * \param domain Random domain. This is used to guarantee that different
594 * applications of a random engine inside the code get different
595 * streams of random numbers, without requiring the user
596 * to provide lots of random seeds. Pick a value from the
597 * RandomDomain class, or RandomDomain::Other if it is
598 * not important. In the latter case you might want to use
599 * \ref gmx::DefaultRandomEngine instead.
601 * \note The random domain is really another 64-bit seed value.
603 * Re-initialized the seed similar to the counter constructor.
604 * Same rules apply: The highest few bits of the last word are
605 * reserved to encode the number of internal counter bits, but
606 * to save the user the trouble of making sure these are zero
607 * when using e.g. a random device, we just ignore them.
609 void seed(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other)
611 seed(key0, static_cast<uint64_t>(domain));
614 /*! \brief Seed random engine from 2x64-bit unsigned integers
616 * This assigns the raw 128 bit key data from unsigned integers.
617 * It is meant for the case when you want full control over the key,
618 * for instance to compare with reference values of the ThreeFry
619 * function during testing.
621 * \param key0 First word of key/random seed.
622 * \param key1 Second word of key/random seed.
624 * \throws InternalError if the high bits needed to encode the number of counter
625 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
627 void seed(uint64_t key0, uint64_t key1)
629 const unsigned int internalCounterBitsBits =
630 (internalCounterBits > 0) ? (StaticLog2<internalCounterBits>::value + 1) : 0;
632 key_ = { { key0, key1 } };
634 if (internalCounterBits > 0)
636 internal::highBitCounter::checkAndClear<result_type, 2, internalCounterBitsBits>(&key_);
637 internal::highBitCounter::increment<result_type, 2, internalCounterBitsBits>(
638 &key_, internalCounterBits - 1);
643 /*! \brief Restart 2x64 random engine counter from 2 64-bit values
645 * \param ctr0 First word of new counter, in the form of 64-bit unsigned values.
646 * \param ctr1 Second word of new counter
648 * Restarting the engine with a new counter is extremely fast with ThreeFry64,
649 * and basically just consists of storing the counter value, so you should
650 * use this liberally in your innermost loops to restart the engine with
651 * e.g. the current step and atom index as counter values.
653 * \throws InternalError if any of the highest bits that are reserved
654 * for the internal part of the counter are set. The number of
655 * reserved bits is to the last template parameter to the class.
657 void restart(uint64_t ctr0 = 0, uint64_t ctr1 = 0)
660 counter_ = { { ctr0, ctr1 } };
661 if (!internal::highBitCounter::checkAndClear<result_type, 2, internalCounterBits>(&counter_))
663 GMX_THROW(InternalError(
664 "High bits of counter are reserved for the internal stream counter."));
666 block_ = generateBlock(key_, counter_);
670 /*! \brief Generate the next random number
672 * This will return the next stored 64-bit value if one is available,
673 * and otherwise generate a new block, update the internal counters, and
674 * return the first value while storing the others.
676 * \throws InternalError if the internal counter space is exhausted.
678 result_type operator()()
680 if (index_ >= c_resultsPerCounter_)
682 internal::highBitCounter::increment<result_type, 2, internalCounterBits>(&counter_);
683 block_ = generateBlock(key_, counter_);
686 return block_[index_++];
689 /*! \brief Skip next n random numbers
691 * Moves the internal random stream for the give key/counter value
692 * n positions forward. The count is based on the number of random values
693 * returned, such that skipping 5 values gives exactly the same result as
694 * drawing 5 values that are ignored.
696 * \param n Number of values to jump forward.
698 * \throws InternalError if the internal counter space is exhausted.
700 void discard(uint64_t n)
702 index_ += n % c_resultsPerCounter_;
703 n /= c_resultsPerCounter_;
705 if (index_ > c_resultsPerCounter_)
707 index_ -= c_resultsPerCounter_;
711 // Make sure the state is the same as if we came to this counter and
712 // index by natural generation.
713 if (index_ == 0 && n > 0)
715 index_ = c_resultsPerCounter_;
718 internal::highBitCounter::increment<result_type, 2, internalCounterBits>(&counter_, n);
719 block_ = generateBlock(key_, counter_);
722 /*! \brief Return true if two ThreeFry2x64 engines are identical
724 * \param x Instance to compare with.
726 * This routine should return true if the two engines will generate
727 * identical random streams when drawing.
729 bool operator==(const ThreeFry2x64General<rounds, internalCounterBits>& x) const
731 // block_ is uniquely specified by key_ and counter_.
732 return (key_ == x.key_ && counter_ == x.counter_ && index_ == x.index_);
735 /*! \brief Return true of two ThreeFry2x64 engines are not identical
737 * \param x Instance to compare with.
739 * This routine should return true if the two engines will generate
740 * different random streams when drawing.
742 bool operator!=(const ThreeFry2x64General<rounds, internalCounterBits>& x) const
744 return !operator==(x);
748 /*! \brief Number of results returned for each invocation of the block generation */
749 static const unsigned int c_resultsPerCounter_ =
750 static_cast<unsigned int>(sizeof(counter_type) / sizeof(result_type));
752 /*! \brief ThreeFry2x64 key, i.e. the random seed for this stream.
754 * The highest few bits of the key are replaced to encode the value of
755 * internalCounterBits, in order to make all streams unique.
759 /*! \brief ThreeFry2x64 total counter.
761 * The highest internalCounterBits are reserved for an internal counter
762 * so that the combination of a key and counter provides a stream that
763 * returns 2*2^internalCounterBits (ThreeFry2x64) random 64-bit values before
764 * the internal counter space is exhausted and an exception is thrown.
766 counter_type counter_;
767 /*! \brief The present block encrypted from values of key and counter. */
769 /*! \brief Index of the next value in block_ to return from random engine */
772 GMX_DISALLOW_COPY_AND_ASSIGN(ThreeFry2x64General);
776 /*! \brief ThreeFry2x64 random engine with 20 iteractions.
778 * \tparam internalCounterBits, default 64.
780 * This class provides very high quality random numbers that pass all
781 * BigCrush tests, it works with two 64-bit values each for keys and
782 * counters, and is most efficient when we only need a few random values
783 * before restarting the counters with new values.
785 template<unsigned int internalCounterBits = 64>
786 class ThreeFry2x64 : public ThreeFry2x64General<20, internalCounterBits>
789 /*! \brief Construct ThreeFry random engine with 2x64 key values, 20 rounds.
791 * \param key0 Random seed in the form of a 64-bit unsigned value.
792 * \param domain Random domain. This is used to guarantee that different
793 * applications of a random engine inside the code get different
794 * streams of random numbers, without requiring the user
795 * to provide lots of random seeds. Pick a value from the
796 * RandomDomain class, or RandomDomain::Other if it is
797 * not important. In the latter case you might want to use
798 * \ref gmx::DefaultRandomEngine instead.
800 * \note The random domain is really another 64-bit seed value.
802 * \throws InternalError if the high bits needed to encode the number of counter
805 ThreeFry2x64(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) :
806 ThreeFry2x64General<20, internalCounterBits>(key0, domain)
810 /*! \brief Construct random engine from 2x64-bit unsigned integers, 20 rounds
812 * This constructor assigns the raw 128 bit key data from unsigned integers.
813 * It is meant for the case when you want full control over the key,
814 * for instance to compare with reference values of the ThreeFry
815 * function during testing.
817 * \param key0 First word of key/random seed.
818 * \param key1 Second word of key/random seed.
820 * \throws InternalError if the high bits needed to encode the number of counter
821 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
823 ThreeFry2x64(uint64_t key0, uint64_t key1) :
824 ThreeFry2x64General<20, internalCounterBits>(key0, key1)
829 /*! \brief ThreeFry2x64 random engine with 13 iteractions.
831 * \tparam internalCounterBits, default 64.
833 * This class provides relatively high quality random numbers that only
834 * fail one BigCrush test, and it is a bit faster than the 20-round version.
835 * It works with two 64-bit values each for keys and counters, and is most
836 * efficient when we only need a few random values before restarting
837 * the counters with new values.
839 template<unsigned int internalCounterBits = 64>
840 class ThreeFry2x64Fast : public ThreeFry2x64General<13, internalCounterBits>
843 /*! \brief Construct ThreeFry random engine with 2x64 key values, 13 rounds.
845 * \param key0 Random seed in the form of a 64-bit unsigned value.
846 * \param domain Random domain. This is used to guarantee that different
847 * applications of a random engine inside the code get different
848 * streams of random numbers, without requiring the user
849 * to provide lots of random seeds. Pick a value from the
850 * RandomDomain class, or RandomDomain::Other if it is
851 * not important. In the latter case you might want to use
852 * \ref gmx::DefaultRandomEngine instead.
854 * \note The random domain is really another 64-bit seed value.
856 * \throws InternalError if the high bits needed to encode the number of counter
859 ThreeFry2x64Fast(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) :
860 ThreeFry2x64General<13, internalCounterBits>(key0, domain)
864 /*! \brief Construct ThreeFry random engine from 2x64-bit unsigned integers, 13 rounds.
866 * This constructor assigns the raw 128 bit key data from unsigned integers.
867 * It is meant for the case when you want full control over the key,
868 * for instance to compare with reference values of the ThreeFry
869 * function during testing.
871 * \param key0 First word of key/random seed.
872 * \param key1 Second word of key/random seed.
874 * \throws InternalError if the high bits needed to encode the number of counter
875 * bits are nonzero. To test arbitrary values, use 0 internal counter bits.
877 ThreeFry2x64Fast(uint64_t key0, uint64_t key1) :
878 ThreeFry2x64General<13, internalCounterBits>(key0, key1)
884 /*! \brief Default fast and accurate random engine in Gromacs
886 * This engine will return 2*2^64 random results using the default
887 * gmx::RandomDomain::Other stream, and can be initialized with a single
888 * seed argument without having to remember empty template angle brackets.
890 typedef ThreeFry2x64Fast<> DefaultRandomEngine;
894 #endif // GMX_RANDOM_THREEFRY_H