libgromacs-dev 2018.1-1 / usr / include / gromacs / random / threefry.h

/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 2015,2016,2017, by the GROMACS development team, led by
 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 * and including many others, as listed in the AUTHORS file in the
 * top-level source directory and at http://www.gromacs.org.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */

/*! \file
 * \brief Implementation of the 2x64 ThreeFry random engine
 *
 * \author Erik Lindahl <erik.lindahl@gmail.com>
 * \inpublicapi
 * \ingroup module_random
 */

#ifndef GMX_RANDOM_THREEFRY_H
#define GMX_RANDOM_THREEFRY_H

#include <array>
#include <limits>

#include "gromacs/math/functions.h"
#include "gromacs/random/seed.h"
#include "gromacs/utility/classhelpers.h"
#include "gromacs/utility/exceptions.h"

/*
 * The GROMACS implementation of the ThreeFry random engine has been
 * heavily inspired by the versions proposed to Boost by:
 *
 * John Salmon, Copyright 2010-2014 by D. E. Shaw Research
 * https://github.com/DEShawResearch/Random123-Boost
 *
 * Thijs van den Berg, Copyright (c) 2014 M.A. (Thijs) van den Berg
 * https://github.com/sitmo/threefry
 *
 * Both of them are covered by the Boost Software License:
 *
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 *
 * The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */

namespace gmx
{

namespace internal
{
// Variable-bitfield counters used to increment internal counters as
// part of std::arrays.

struct
highBitCounter
{
    /*! \brief Clear highBits higest bits of ctr, return false if they were non-zero.
     *
     *  This function clears the space required for the internal counters,
     *  and returns true if they were correctly zero when calling, false otherwise.
     *
     *  \tparam        UIntType  Integer type to use for each word in counter
     *  \tparam        words     Number of UIntType words in counter
     *  \tparam        highBits  Number of bits to check. The template parameter makes it
     *                           possible to optimize this extensively at compile time.
     *  \param         ctr       Reference to counter to check and clear.
     */
    template<class UIntType, std::size_t words, unsigned int highBits>
    static bool
    checkAndClear(std::array<UIntType, words> * ctr)
    {
        const std::size_t  bitsPerWord       = std::numeric_limits<UIntType>::digits;
        const std::size_t  bitsTotal         = bitsPerWord*words;

        static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");

        const std::size_t  lastWordIdx       = (bitsTotal - highBits) / bitsPerWord;
        const std::size_t  lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
        const UIntType     lastWordOne       = static_cast<UIntType>(1) << lastWordLowBitIdx;
        const UIntType     mask              = lastWordOne-1;

        bool               isClear              = true;

        for (unsigned int i = words-1; i > lastWordIdx; --i)
        {
            if ((*ctr)[i])
            {
                isClear    = false;
                (*ctr)[i]  = 0;
            }
        }
        if (highBits > 0 && (*ctr)[lastWordIdx] >= lastWordOne)
        {
            isClear                 = false;
            (*ctr)[lastWordIdx]    &= mask;
        }
        return isClear;
    }

    /*! \brief Increment the internal counter in highBits by one
     *
     *  \tparam         UIntType  Integer type to use for each word in counter
     *  \tparam         words     Number of UIntType words in counter
     *  \tparam         highBits  Number of bits reserved for the internal counter.
     *  \param          ctr       Reference to the counter value to increment.
     *
     *  \throws InternalError if internal counter space is exhausted.
     *
     *  This routine will work across the word boundaries for any number
     *  of internal counter bits that fits in the total counter.
     */
    template<class UIntType, std::size_t words, unsigned int highBits>
    static void
    increment(std::array<UIntType, words> * ctr)
    {
        const std::size_t  bitsPerWord       = std::numeric_limits<UIntType>::digits;
        const std::size_t  bitsTotal         = bitsPerWord*words;

        static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");

        const std::size_t  lastWordIdx       = (bitsTotal - highBits) / bitsPerWord;
        const std::size_t  lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
        const UIntType     lastWordOne       = static_cast<UIntType>(1) << lastWordLowBitIdx;

        // For algorithm & efficiency reasons we need to store the internal counter in
        // the same array as the user-provided counter, so we use the higest bits, possibly
        // crossing several words.
        //
        // To have the computer help us with the dirty carry arithmetics we store the bits
        // in the internal counter part in normal fashion, but the internal counter words in
        // reverse order; the highest word of the total counter array (words-1) is thus
        // the least significant part of the internal counter (if it spans several words).
        //
        // The incrementation works as follows:
        //
        // 0) If the index of the least significant internal counter word is larger
        //    than words-1, there was never any space.
        // 1) If the internal counter spans more than one word, we must have one or
        //    more internal counter words that correspond entirely to the this counter.
        //    Start with the least significant one (words-1) and increment it.
        //    If the new value is not zero we did not loop around (no carry), so everything
        //    is good, and we are done - return!
        //    If the new value is zero, we need to move the carry result to the next word,
        //    so we just continue the loop until we have gone through all words that
        //    are internal-counter-only.
        // 2) After the loop, there is stuff remaining to add, and by definition there
        //    is some internal counter space in the next word, but the question
        //    is if we have exhausted it. We already created a constant that corresponds
        //    to the bit that represents '1' for the internal counter part of this word.
        //    When we add this constant it will not affect the user-counter-part at all,
        //    and if we exhaust the internal counter space the high bits will cause the entire
        //    word to wrap around, and the result will be smaller than the bit we added.
        //    If this happens we throw, otherwise we're done.
        //
        // Since all constants will be evaluated at compile time, this entire routine
        // will usually be reduced to simply incrementing a word by a constant, and throwing
        // if the result is smaller than the constant.

        if (lastWordIdx >= words)
        {
            GMX_THROW(InternalError("Cannot increment random engine defined with 0 internal counter bits."));
        }

        for (unsigned int i = words-1; i > lastWordIdx; --i)
        {
            (*ctr)[i]++;
            if ((*ctr)[i])
            {
                return;     // No carry means we are done
            }
        }
        (*ctr)[lastWordIdx] += lastWordOne;
        if ((*ctr)[lastWordIdx] < lastWordOne)
        {
            GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
        }
        return;
    }

    /*! \brief Increment the internal counter in highBits by a value.
     *
     *  \tparam        UIntType  Integer type to use for each word in counter
     *  \tparam        words     Number of UIntType words in counter
     *  \tparam        highBits  Number of bits reserved for the internal counter.
     *  \param         ctr       Reference to the counter to increment.
     *  \param         addend    Value to add to internal.
     *
     *  \throws InternalError if internal counter space is exhausted.
     *
     *  This routine will work across the word boundaries for any number
     *  of internal counter bits that fits in the total counter.
     */
    template<class UIntType, std::size_t words, unsigned int highBits>
    static void
    increment(std::array<UIntType, words> * ctr, UIntType addend)
    {
        const std::size_t  bitsPerWord       = std::numeric_limits<UIntType>::digits;
        const std::size_t  bitsTotal         = bitsPerWord*words;

        static_assert(highBits <= bitsTotal, "High bits do not fit in counter.");

        const std::size_t  lastWordIdx       = (bitsTotal - highBits) / bitsPerWord;
        const std::size_t  lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord;
        const UIntType     lastWordOne       = static_cast<UIntType>(1) << lastWordLowBitIdx;
        const UIntType     lastWordMaxVal    = (~static_cast<UIntType>(0)) >> lastWordLowBitIdx;

        if (lastWordIdx >= words)
        {
            GMX_THROW(InternalError("Cannot increment random engine defined with 0 internal counter bits."));
        }

        for (unsigned int i = words-1; i > lastWordIdx; --i)
        {
            (*ctr)[i] += addend;
            addend     = ((*ctr)[i] < addend);   // 1 is the carry!
            if (addend == 0)
            {
                return;
            }
        }

        if (addend > lastWordMaxVal)
        {
            GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
        }
        addend *= lastWordOne;

        (*ctr)[lastWordIdx] += addend;

        if ((*ctr)[lastWordIdx] < addend)
        {
            GMX_THROW(InternalError("Random engine stream ran out of internal counter space."));
        }
        return;
    }
};
}

/*! \brief General implementation class for ThreeFry counter-based random engines.
 *
 *  This class is used to implement several different ThreeFry2x64 random engines
 *  differing in the number of rounds executed in and the number of bits reserved
 *  for the internal counter. It is compatible with C++11 random engines, and
 *  can be used e.g. with all random distributions from the standard library.
 *
 *  ThreeFry is a counter-based rather than state-based random engine. This
 *  means that we seed it with a "key", after which we can get the
 *  N:th random number in a sequence (specified by a counter) directly. This
 *  means we are guaranteed the same sequence of numbers even when running in
 *  parallel if using e.g. step and atom index as counters.
 *
 *  However, it is also useful to be able to use it as a normal random engine,
 *  for instance if you need more than 2 64-bit random values for a specific
 *  counter value, not to mention where you just need good normal random numbers.
 *  To achieve this, this implementation uses John Salmon's idea of reserving
 *  a couple of the highest bits in the user-provided counter for an internal
 *  counter. For instance, if reserving 3 bits, this means you get a stream of
 *  8 iterations (each with 2 random values) after every restart. If you call
 *  the engine after these bits have been exhausted, it will throw an
 *  exception to make sure you don't get overlapping streams by mistake.
 *  Reserving 3 bits also means you can only use 64-3=61 bits of the highest
 *  word when restarting (i.e., setting) the counters.
 *
 *  This version also supports using internalCounterBits=0. In this case the
 *  random engine will be able to return a single counter round, i.e. 2 64-bit
 *  values for ThreeFry2x64, after which an exception is thrown. In this case no
 *  high bits are reserved, which means the class implements the raw ThreeFry2x64
 *  random function.
 *
 *  \tparam rounds  The number of encryption iterations used when generating.
 *                  This can in principle be any value, but 20 rounds has been
 *                  shown to pass all BigCrush random tests, and with 13 rounds
 *                  only one fails. This is a very stringent test, and the
 *                  standard Mersenne Twister engine fails two, so 13 rounds
 *                  should be a perfectly fine balance in most cases.
 *  \tparam internalCounterBits
 *                  Number of high bits in the user-provided counter reserved
 *                  for the internal counter. The number of values the engine
 *                  can return after each restart will be
 *                  words*2^internalCounterBits.
 */
template<unsigned int rounds, unsigned int internalCounterBits>
class ThreeFry2x64General
{
    // While this class will formally work with any value for rounds, there is
    // no reason to go lower than 13, and this might help catch some typos.
    // If we find a reason to use lower values in the future, or if you simply
    // want to test, this assert can safely be removed.
    static_assert(rounds >= 13, "You should not use less than 13 encryption rounds for ThreeFry2x64.");

    public:
        // result_type must be lower case to be compatible with C++11 standard library

        /*! \brief Integer type for output. */
        typedef gmx_uint64_t                    result_type;
        /*! \brief Use array for counter & key states so it is allocated on the stack */
        typedef std::array<result_type, 2>      counter_type;

    private:

        /*! \brief Rotate value left by specified number of bits
         *
         *  \param i    Value to rotate (result_type, which should be 64-bit).
         *  \param bits Number of bits to rotate i.
         *
         *  \return Input value rotated 'bits' left.
         */
        result_type
        rotLeft(result_type i, unsigned int bits)
        {
            return (i << bits) | (i >> (std::numeric_limits<result_type>::digits-bits));
        }

        /*! \brief Perform encryption step for ThreeFry2x64 algorithm
         *
         *  It performs the encryption step of the standard ThreeFish symmetric-key
         *  tweakable block cipher, which is the core of the ThreeFry random
         *  engine. The number of encryption rounds is specified by the class
         *  template parameter 'rounds'.
         *
         *  \param key   Reference to key value
         *  \param ctr   Counter value to use
         *
         *  \return Newly encrypted 2x64 block, according to the class template parameters.
         */
        counter_type
        generateBlock(const counter_type &key,
                      const counter_type &ctr)
        {
            const unsigned int  rotations[] = {16, 42, 12, 31, 16, 32, 24, 21};
            counter_type        x           = ctr;

            result_type         ks[3] = { 0x0, 0x0, 0x1bd11bdaa9fc1a22 };

            // This is actually a pretty simple routine that merely executes the
            // for-block specified further down 'rounds' times. However, both
            // clang and gcc have problems unrolling and replacing rotations[r%8]
            // with constants, so we unroll the first 20 iterations manually.

            if (rounds > 0)
            {
                ks[0] = key[0]; ks[2] ^= key[0]; x[0] = x[0] + key[0];
                ks[1] = key[1]; ks[2] ^= key[1]; x[1] = x[1] + key[1];
                x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0];
            }
            if (rounds > 1)  { x[0] += x[1]; x[1] = rotLeft(x[1], 42); x[1] ^= x[0]; }
            if (rounds > 2)  { x[0] += x[1]; x[1] = rotLeft(x[1], 12); x[1] ^= x[0]; }
            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; }
            if (rounds > 4)  { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
            if (rounds > 5)  { x[0] += x[1]; x[1] = rotLeft(x[1], 32); x[1] ^= x[0]; }
            if (rounds > 6)  { x[0] += x[1]; x[1] = rotLeft(x[1], 24); x[1] ^= x[0]; }
            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; }
            if (rounds > 8)  { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
            if (rounds > 9)  { x[0] += x[1]; x[1] = rotLeft(x[1], 42); x[1] ^= x[0]; }
            if (rounds > 10) { x[0] += x[1]; x[1] = rotLeft(x[1], 12); x[1] ^= x[0]; }
            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; }
            if (rounds > 12) { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
            if (rounds > 13) { x[0] += x[1]; x[1] = rotLeft(x[1], 32); x[1] ^= x[0]; }
            if (rounds > 14) { x[0] += x[1]; x[1] = rotLeft(x[1], 24); x[1] ^= x[0]; }
            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; }
            if (rounds > 16) { x[0] += x[1]; x[1] = rotLeft(x[1], 16); x[1] ^= x[0]; }
            if (rounds > 17) { x[0] += x[1]; x[1] = rotLeft(x[1], 42); x[1] ^= x[0]; }
            if (rounds > 18) { x[0] += x[1]; x[1] = rotLeft(x[1], 12); x[1] ^= x[0]; }
            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; }

            for (unsigned int r = 20; r < rounds; r++)
            {
                x[0] += x[1];
                x[1]  = rotLeft(x[1], rotations[r%8]);
                x[1] ^= x[0];
                if (( (r + 1) & 3 ) == 0)
                {
                    unsigned int r4 = (r + 1) >> 2;
                    x[0] += ks[ r4 % 3 ];
                    x[1] += ks[ (r4 + 1) % 3 ] + r4;
                }
            }
            return x;
        }

    public:
        //! \brief Smallest value that can be returned from random engine.
#if !defined(_MSC_VER)
        static constexpr
#else
        // Avoid constexpr bug in MSVC 2015, note that max() below does work
        static
#endif
        result_type min() { return std::numeric_limits<result_type>::min(); }

        //! \brief Largest value that can be returned from random engine.
        static constexpr
        result_type max() { return std::numeric_limits<result_type>::max(); }

        /*! \brief Construct random engine with 2x64 key values
         *
         *  This constructor takes two values, and should only be used with
         *  the 2x64 implementations.
         *
         *  \param key0   Random seed in the form of a 64-bit unsigned value.
         *  \param domain Random domain. This is used to guarantee that different
         *                applications of a random engine inside the code get different
         *                streams of random numbers, without requiring the user
         *                to provide lots of random seeds. Pick a value from the
         *                RandomDomain class, or RandomDomain::Other if it is
         *                not important. In the latter case you might want to use
         *                \ref gmx::DefaultRandomEngine instead.
         *
         *  \note The random domain is really another 64-bit seed value.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero.
         */
        ThreeFry2x64General(gmx_uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other)
        {
            seed(key0, domain);
        }

        /*! \brief Construct random engine from 2x64-bit unsigned integers
         *
         *  This constructor assigns the raw 128 bit key data from unsigned integers.
         *  It is meant for the case when you want full control over the key,
         *  for instance to compare with reference values of the ThreeFry
         *  function during testing.
         *
         *  \param key0   First word of key/random seed.
         *  \param key1   Second word of key/random seed.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero. To test arbitrary values, use 0 internal counter bits.
         */
        ThreeFry2x64General(gmx_uint64_t key0, gmx_uint64_t key1)
        {
            seed(key0, key1);
        }

        /*! \brief Seed 2x64 random engine with two 64-bit key values
         *
         *  \param key0   First word of random seed, in the form of 64-bit unsigned values.
         *  \param domain Random domain. This is used to guarantee that different
         *                applications of a random engine inside the code get different
         *                streams of random numbers, without requiring the user
         *                to provide lots of random seeds. Pick a value from the
         *                RandomDomain class, or RandomDomain::Other if it is
         *                not important. In the latter case you might want to use
         *                \ref gmx::DefaultRandomEngine instead.
         *
         *  \note The random domain is really another 64-bit seed value.
         *
         *  Re-initialized the seed similar to the counter constructor.
         *  Same rules apply: The highest few bits of the last word are
         *  reserved to encode the number of internal counter bits, but
         *  to save the user the trouble of making sure these are zero
         *  when using e.g. a random device, we just ignore them.
         */
        void
        seed(gmx_uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other)
        {
            seed(key0, static_cast<gmx_uint64_t>(domain));
        }

        /*! \brief Seed random engine from 2x64-bit unsigned integers
         *
         *  This assigns the raw 128 bit key data from unsigned integers.
         *  It is meant for the case when you want full control over the key,
         *  for instance to compare with reference values of the ThreeFry
         *  function during testing.
         *
         *  \param key0   First word of key/random seed.
         *  \param key1   Second word of key/random seed.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero. To test arbitrary values, use 0 internal counter bits.
         */
        void
        seed(gmx_uint64_t key0, gmx_uint64_t key1)
        {
            const unsigned int internalCounterBitsBits = (internalCounterBits > 0) ? ( StaticLog2<internalCounterBits>::value + 1 ) : 0;

            key_ = {{key0, key1}};

            if (internalCounterBits > 0)
            {
                internal::highBitCounter::checkAndClear<result_type, 2, internalCounterBitsBits>(&key_);
                internal::highBitCounter::increment<result_type, 2, internalCounterBitsBits>(&key_, internalCounterBits-1);
            }
            restart(0, 0);
        }

        /*! \brief Restart 2x64 random engine counter from 2 64-bit values
         *
         *  \param ctr0 First word of new counter, in the form of 64-bit unsigned values.
         *  \param ctr1 Second word of new counter
         *
         * Restarting the engine with a new counter is extremely fast with ThreeFry64,
         * and basically just consists of storing the counter value, so you should
         * use this liberally in your innermost loops to restart the engine with
         * e.g. the current step and atom index as counter values.
         *
         * \throws InternalError if any of the highest bits that are reserved
         *         for the internal part of the counter are set. The number of
         *         reserved bits is to the last template parameter to the class.
         */
        void
        restart(gmx_uint64_t ctr0 = 0, gmx_uint64_t ctr1 = 0)
        {

            counter_ = {{ctr0, ctr1}};
            if (!internal::highBitCounter::checkAndClear<result_type, 2, internalCounterBits>(&counter_))
            {
                GMX_THROW(InternalError("High bits of counter are reserved for the internal stream counter."));
            }
            block_ = generateBlock(key_, counter_);
            index_ = 0;
        }

        /*! \brief Generate the next random number
         *
         *  This will return the next stored 64-bit value if one is available,
         *  and otherwise generate a new block, update the internal counters, and
         *  return the first value while storing the others.
         *
         *  \throws InternalError if the internal counter space is exhausted.
         */
        result_type
        operator()()
        {
            if (index_ >= c_resultsPerCounter_)
            {
                internal::highBitCounter::increment<result_type, 2, internalCounterBits>(&counter_);
                block_ = generateBlock(key_, counter_);
                index_ = 0;
            }
            return block_[index_++];
        }

        /*! \brief Skip next n random numbers
         *
         *  Moves the internal random stream for the give key/counter value
         *  n positions forward. The count is based on the number of random values
         *  returned, such that skipping 5 values gives exactly the same result as
         *  drawing 5 values that are ignored.
         *
         *  \param n Number of values to jump forward.
         *
         *  \throws InternalError if the internal counter space is exhausted.
         */
        void
        discard(gmx_uint64_t n)
        {
            index_ += n % c_resultsPerCounter_;
            n      /= c_resultsPerCounter_;

            if (index_ > c_resultsPerCounter_)
            {
                index_ -= c_resultsPerCounter_;
                n++;
            }

            // Make sure the state is the same as if we came to this counter and
            // index by natural generation.
            if (index_ == 0 && n > 0)
            {
                index_ = c_resultsPerCounter_;
                n--;
            }
            internal::highBitCounter::increment<result_type, 2, internalCounterBits>(&counter_, n);
            block_ = generateBlock(key_, counter_);
        }

        /*! \brief Return true if two ThreeFry2x64 engines are identical
         *
         * \param  x    Instance to compare with.
         *
         * This routine should return true if the two engines will generate
         * identical random streams when drawing.
         */
        bool
        operator==(const ThreeFry2x64General<rounds, internalCounterBits> &x) const
        {
            // block_ is uniquely specified by key_ and counter_.
            return (key_ == x.key_ && counter_ == x.counter_ && index_ == x.index_);
        }

        /*! \brief Return true of two ThreeFry2x64 engines are not identical
         *
         * \param  x    Instance to compare with.
         *
         * This routine should return true if the two engines will generate
         * different random streams when drawing.
         */
        bool
        operator!=(const ThreeFry2x64General<rounds, internalCounterBits> &x) const { return !operator==(x); }

    private:

        /*! \brief Number of results returned for each invocation of the block generation */
        static const unsigned int c_resultsPerCounter_  = static_cast<unsigned int>(sizeof(counter_type)/sizeof(result_type));

        /*! \brief ThreeFry2x64 key, i.e. the random seed for this stream.
         *
         *  The highest few bits of the key are replaced to encode the value of
         *  internalCounterBits, in order to make all streams unique.
         */
        counter_type key_;

        /*! \brief ThreeFry2x64 total counter.
         *
         *  The highest internalCounterBits are reserved for an internal counter
         *  so that the combination of a key and counter provides a stream that
         *  returns 2*2^internalCounterBits (ThreeFry2x64) random 64-bit values before
         *  the internal counter space is exhausted and an exception is thrown.
         */
        counter_type counter_;
        /*! \brief The present block encrypted from values of key and counter. */
        counter_type block_;
        /*! \brief Index of the next value in block_ to return from random engine */
        unsigned int index_;

        GMX_DISALLOW_COPY_AND_ASSIGN(ThreeFry2x64General);
};


/*! \brief ThreeFry2x64 random engine with 20 iteractions.
 *
 *  \tparam internalCounterBits, default 64.
 *
 *  This class provides very high quality random numbers that pass all
 *  BigCrush tests, it works with two 64-bit values each for keys and
 *  counters, and is most  efficient when we only need a few random values
 *  before restarting the counters with new values.
 */
template<unsigned int internalCounterBits = 64>
class ThreeFry2x64 : public ThreeFry2x64General<20, internalCounterBits>
{
    public:
        /*! \brief Construct ThreeFry random engine with 2x64 key values, 20 rounds.
         *
         *  \param key0   Random seed in the form of a 64-bit unsigned value.
         *  \param domain Random domain. This is used to guarantee that different
         *                applications of a random engine inside the code get different
         *                streams of random numbers, without requiring the user
         *                to provide lots of random seeds. Pick a value from the
         *                RandomDomain class, or RandomDomain::Other if it is
         *                not important. In the latter case you might want to use
         *                \ref gmx::DefaultRandomEngine instead.
         *
         *  \note The random domain is really another 64-bit seed value.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero.
         */
        ThreeFry2x64(gmx_uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) : ThreeFry2x64General<20, internalCounterBits>(key0, domain) {}

        /*! \brief Construct random engine from 2x64-bit unsigned integers, 20 rounds
         *
         *  This constructor assigns the raw 128 bit key data from unsigned integers.
         *  It is meant for the case when you want full control over the key,
         *  for instance to compare with reference values of the ThreeFry
         *  function during testing.
         *
         *  \param key0   First word of key/random seed.
         *  \param key1   Second word of key/random seed.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero. To test arbitrary values, use 0 internal counter bits.
         */
        ThreeFry2x64(gmx_uint64_t key0, gmx_uint64_t key1) : ThreeFry2x64General<20, internalCounterBits>(key0, key1) {}
};

/*! \brief ThreeFry2x64 random engine with 13 iteractions.
 *
 *  \tparam internalCounterBits, default 64.
 *
 *  This class provides relatively high quality random numbers that only
 *  fail one BigCrush test, and it is a bit faster than the 20-round version.
 *  It works with two 64-bit values each for keys and counters, and is most
 *  efficient when we only need a few random values before restarting
 *  the counters with new values.
 */
template<unsigned int internalCounterBits = 64>
class ThreeFry2x64Fast : public ThreeFry2x64General<13, internalCounterBits>
{
    public:
        /*! \brief Construct ThreeFry random engine with 2x64 key values, 13 rounds.
         *
         *  \param key0   Random seed in the form of a 64-bit unsigned value.
         *  \param domain Random domain. This is used to guarantee that different
         *                applications of a random engine inside the code get different
         *                streams of random numbers, without requiring the user
         *                to provide lots of random seeds. Pick a value from the
         *                RandomDomain class, or RandomDomain::Other if it is
         *                not important. In the latter case you might want to use
         *                \ref gmx::DefaultRandomEngine instead.
         *
         *  \note The random domain is really another 64-bit seed value.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero.
         */
        ThreeFry2x64Fast(gmx_uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) : ThreeFry2x64General<13, internalCounterBits>(key0, domain) {}

        /*! \brief Construct ThreeFry random engine from 2x64-bit unsigned integers, 13 rounds.
         *
         *  This constructor assigns the raw 128 bit key data from unsigned integers.
         *  It is meant for the case when you want full control over the key,
         *  for instance to compare with reference values of the ThreeFry
         *  function during testing.
         *
         *  \param key0   First word of key/random seed.
         *  \param key1   Second word of key/random seed.
         *
         *  \throws InternalError if the high bits needed to encode the number of counter
         *          bits are nonzero. To test arbitrary values, use 0 internal counter bits.
         */
        ThreeFry2x64Fast(gmx_uint64_t key0, gmx_uint64_t key1) : ThreeFry2x64General<13, internalCounterBits>(key0, key1) {}
};



/*! \brief Default fast and accurate random engine in Gromacs
 *
 *  This engine will return 2*2^64 random results using the default
 *  gmx::RandomDomain::Other stream, and can be initialized with a single
 *  seed argument without having to remember empty template angle brackets.
 */
typedef ThreeFry2x64Fast<>                  DefaultRandomEngine;

}      // namespace gmx

#endif // GMX_RANDOM_THREEFRY_H