/usr/include/gromacs/random/threefry.h is in libgromacs-dev 2016.1-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 | /*
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2015,2016, 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.
static gmx_constexpr
result_type min() { return std::numeric_limits<result_type>::min(); }
//! \brief Largest value that can be returned from random engine.
static gmx_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
|