diff --git a/README.md b/README.md index b08c09b2a2..b8894d233b 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ Fast, parallelized molecular dynamics trajectory data analysis. Build Status ============= -* Travis-CI: [![Travis Build Status](https://travis-ci.org/Amber-MD/cpptraj.svg?branch=master)](https://travis-ci.org/Amber-MD/cpptraj) +* GitHub Actions: [![GitHub Actions Status](https://github.com/Amber-MD/cpptraj/actions/workflows/merge-gate.yml/badge.svg)](https://github.com/Amber-MD/cpptraj/actions) * AppVeyor: [![Windows Build Status](https://ci.appveyor.com/api/projects/status/github/Amber-MD/cpptraj?branch=master&svg=true&retina=true)](https://ci.appveyor.com/project/drroe/cpptraj-aof9y/branch/master) * Jenkins: [![Jenkins Build Status](https://jenkins.jasonswails.com/buildStatus/icon?job=amber-github%2Fcpptraj%2Fmaster&style=plastic)](https://jenkins.jasonswails.com/job/amber-github/job/cpptraj/job/master/) * LGTM: [![Language grade: C/C++](https://img.shields.io/lgtm/grade/cpp/g/Amber-MD/cpptraj.svg?logo=lgtm&logoWidth=18)](https://lgtm.com/projects/g/Amber-MD/cpptraj/context:cpp) @@ -84,10 +84,12 @@ the following libraries: CPPTRAJ also makes use of the following libraries that are bundled with CPPTRAJ. External ones can be used in place of these if desired. -* ARPACK; without this diagonalization of sparse matrices in `diagmatrix` will be slow -* [helPME](https://github.com/andysim/helpme), required for PME functionality -* XDR for reading GROMACS XTC trajectories -* TNG for reading GROMACS TNG trajectories +* ARPACK; without this diagonalization of sparse matrices in `diagmatrix` will be slow. +* [helPME](https://github.com/andysim/helpme) by Andy Simmonett, required for PME functionality. +* XDR for reading GROMACS XTC trajectories. +* TNG for reading GROMACS TNG trajectories. + +CPPTRAJ also uses the PCG32 and Xoshiro 128++ pseudo-random number generators. `./configure gnu` should be adequate to set up compilation for most systems. For systems without BLAS/LAPACK/ARPACK and/or NetCDF libraries installed, @@ -213,8 +215,8 @@ Original implementation of the Amber NetCDF trajectory format. * Hannes H. Loeffler (STFC Daresbury, Scientific Computing Department, Warrington, WA4 4AD, UK) Diffusion calculation code adapted for use in Action\_STFC\_Diffusion. -External libraries bundled with CPPTRAJ -======================================= +External code/libraries bundled with CPPTRAJ +============================================ * CPPTRAJ makes use of the [GNU readline](https://tiswww.case.edu/php/chet/readline/rltop.html) library for the interactive command line. @@ -227,3 +229,5 @@ External libraries bundled with CPPTRAJ * The reciprocal part of the PME calculation is handled by the [helPME](https://github.com/andysim/helpme) library by Andy Simmonett. * Support for reading DTR trajectories uses the VMD DTR plugin. + +* CPPTRAJ uses code for the [permuted congruent pseudo-random number generator](https://www.pcg-random.org/index.html) PCG32 by Melissa O'Neill and the [Xoshiro 128++ pseudo-random number generator](http://prng.di.unimi.it) by David Blackman and Sebastino Vigna. diff --git a/devtools/ci/jenkins/Dockerfile.libcpptraj b/devtools/ci/jenkins/Dockerfile.libcpptraj index 7cc161264d..10c561fcbd 100644 --- a/devtools/ci/jenkins/Dockerfile.libcpptraj +++ b/devtools/ci/jenkins/Dockerfile.libcpptraj @@ -1,10 +1,11 @@ FROM condaforge/miniforge3 +# noninteractive to try and fix tzdata issue RUN mkdir /cpptraj && \ mkdir /.conda /.local && \ apt-get update && \ apt-get upgrade -y && \ - apt-get install -y gfortran gcc g++ bzip2 libfftw3-dev automake make libbz2-dev \ + DEBIAN_FRONTEND="noninteractive" apt-get install -y gfortran gcc g++ bzip2 libfftw3-dev automake make libbz2-dev \ mpich libmpich-dev zlib1g-dev netcdf-bin libnetcdf-dev \ liblapack-dev libarpack2-dev diff --git a/devtools/rngtest/Makefile b/devtools/rngtest/Makefile new file mode 100644 index 0000000000..40797c11bd --- /dev/null +++ b/devtools/rngtest/Makefile @@ -0,0 +1,47 @@ +include ../../config.h + +#OBJECTS=rngtest.o Random.o RNG.o RNG_Stdlib.o RNG_Marsaglia.o CpptrajStdio.o RNG_MersenneTwister.o RNG_PCG32.o xoshiro128plusplus.o RNG_Xoshiro128pp.o +OBJECTS=xoshiro128plusplus.o +SOURCES=rngtest.cpp ../../src/RNG_MersenneTwister.cpp ../../src/RNG.cpp ../../src/CpptrajStdio.cpp ../../src/Random.cpp ../../src/RNG_Stdlib.cpp ../../src/RNG_Marsaglia.cpp ../../src/RNG_PCG32.cpp ../../src/RNG_Xoshiro128pp.cpp + +all: rngtest + +debug: + make CXXFLAGS='-O0 -g' CFLAGS='-O0 -g' rngtest + +rngtest: $(SOURCES) $(OBJECTS) + $(CXX) $(DIRECTIVES) $(CXXFLAGS) -o rngtest $(SOURCES) $(OBJECTS) + +clean: + /bin/rm -f $(OBJECTS) rngtest + +rngtest.o: rngtest.cpp + $(CXX) $(CXXFLAGS) -c -o rngtest.o rngtest.cpp + +Random.o: ../../src/Random.cpp + $(CXX) $(CXXFLAGS) -c -o Random.o ../../src/Random.cpp + +RNG.o: ../../src/RNG.cpp + $(CXX) $(CXXFLAGS) -c -o RNG.o ../../src/RNG.cpp + +RNG_Stdlib.o: ../../src/RNG_Stdlib.cpp + $(CXX) $(CXXFLAGS) -c -o RNG_Stdlib.o ../../src/RNG_Stdlib.cpp + +RNG_Marsaglia.o: ../../src/RNG_Marsaglia.cpp + $(CXX) $(CXXFLAGS) -c -o RNG_Marsaglia.o ../../src/RNG_Marsaglia.cpp + +RNG_MersenneTwister.o: ../../src/RNG_MersenneTwister.cpp + $(CXX) $(DIRECTIVES) $(CXXFLAGS) -c -o RNG_MersenneTwister.o ../../src/RNG_MersenneTwister.cpp + +RNG_PCG32.o: ../../src/RNG_PCG32.cpp + $(CXX) $(CXXFLAGS) -c -o RNG_PCG32.o ../../src/RNG_PCG32.cpp + +CpptrajStdio.o: ../../src/CpptrajStdio.cpp + $(CXX) $(CXXFLAGS) -c -o CpptrajStdio.o ../../src/CpptrajStdio.cpp + +RNG_Xoshiro128pp.o: ../../src/RNG_Xoshiro128pp.cpp + $(CXX) $(CXXFLAGS) -c -o RNG_Xoshiro128pp.o ../../src/RNG_Xoshiro128pp.cpp + +xoshiro128plusplus.o: ../../src/xoshiro128plusplus.c + $(CC) $(CFLAGS) -c -o xoshiro128plusplus.o ../../src/xoshiro128plusplus.c + diff --git a/devtools/rngtest/README.md b/devtools/rngtest/README.md new file mode 100644 index 0000000000..ea80334d36 --- /dev/null +++ b/devtools/rngtest/README.md @@ -0,0 +1,41 @@ +Dieharder Results +================= + +The results from the Dieharder test suite for the various RNGs in cpptraj are in the results.X files. The following tests were run: +``` + Test Number Test Name Test Reliability +=============================================================================== + -d 0 Diehard Birthdays Test Good + -d 1 Diehard OPERM5 Test Good + -d 2 Diehard 32x32 Binary Rank Test Good + -d 3 Diehard 6x8 Binary Rank Test Good + -d 4 Diehard Bitstream Test Good + -d 8 Diehard Count the 1s (stream) Test Good + -d 9 Diehard Count the 1s Test (byte) Good + -d 10 Diehard Parking Lot Test Good + -d 11 Diehard Minimum Distance (2d Circle) Test Good + -d 12 Diehard 3d Sphere (Minimum Distance) Test Good + -d 13 Diehard Squeeze Test Good + -d 15 Diehard Runs Test Good + -d 16 Diehard Craps Test Good + -d 17 Marsaglia and Tsang GCD Test Good + -d 100 STS Monobit Test Good + -d 101 STS Runs Test Good + -d 102 STS Serial Test (Generalized) Good + -d 201 RGB Generalized Minimum Distance Test Good + -d 202 RGB Permutations Test Good + -d 203 RGB Lagged Sum Test Good + -d 204 RGB Kolmogorov-Smirnov Test Test Good + -d 205 Byte Distribution Good + -d 206 DAB DCT Good + -d 207 DAB Fill Tree Test Good + -d 208 DAB Fill Tree 2 Test Good + -d 209 DAB Monobit 2 Test Good +``` + +- 0: Marsaglia: 44 failed +- 1: C stdlib: 37 failed +- 2: Mersenne Twister (mt19937 from C++ std): 2 failed. +- 3: Permuted congruential generator 32: 2 failed. +- 4: Xoshiro128++: 2 failed. + diff --git a/devtools/rngtest/RunTest.sh b/devtools/rngtest/RunTest.sh new file mode 100755 index 0000000000..00b22b15df --- /dev/null +++ b/devtools/rngtest/RunTest.sh @@ -0,0 +1,50 @@ +#!/bin/bash + +COUNT=100000000 + +make rngtest + +# NOTE: Omitting test 200 for now. Gives following message: +#Error: Can only test distribution of positive ntuples. +# Use -n ntuple for 0 < ntuple. +# Read test description with dieharder -d 200 -h. + +TESTNUMS='0 1 2 3 4 8 9 10 11 12 13 15 16 17 100 101 102 201 202 203 204 205 206 207 208 209' + +TYPELIST='' +DIEHARD='no' +while [ ! -z "$1" ] ; do + case "$1" in + '-type' ) shift; TYPELIST="$TYPELIST $1" ;; + esac + shift +done +if [ -z "$TYPELIST" ] ; then + TYPELIST='0 1 2 3 4' + DIEHARD='yes' +fi +echo "Type list: $TYPELIST" + +for CPPTRAJTYPE in $TYPELIST ; do + OUTFILE=results.$CPPTRAJTYPE + echo "Cpptraj type $CPPTRAJTYPE" + echo "Cpptraj type $CPPTRAJTYPE" > $OUTFILE + TMPFILE=temp.$CPPTRAJTYPE + ./rngtest -t $COUNT -S 1 -r $CPPTRAJTYPE > $TMPFILE + for TEST in $TESTNUMS ; do + echo "Test $TEST" + dieharder -d $TEST -f $TMPFILE -g 202 >> $OUTFILE + done + echo "" +done + +if [ "$DIEHARD" = 'yes' ] ; then + echo "Diehard numbers" + echo "Diehard numbers" > results.diehard + dieharder -S 1 -B -o -t $COUNT > diehard.numbers + for TEST in $TESTNUMS ; do + echo "Test $TEST" + dieharder -d $TEST -f diehard.numbers -g 202 >> results.diehard + done + echo "" +fi diff --git a/devtools/rngtest/results.0 b/devtools/rngtest/results.0 new file mode 100644 index 0000000000..056d2c51a1 --- /dev/null +++ b/devtools/rngtest/results.0 @@ -0,0 +1,260 @@ +Cpptraj type 0 +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.94e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_birthdays| 0| 100| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.08e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_operm5| 0| 1000000| 100|0.74056655| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.81e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_32x32| 0| 40000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.74e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_6x8| 0| 100000| 100|0.94608121| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.69e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_str| 0| 256000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.81e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_byt| 0| 256000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.10e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_parking_lot| 0| 12000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.11e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_2dsphere| 2| 8000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.10e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_3dsphere| 3| 4000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.12e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_squeeze| 0| 100000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.13e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_runs| 0| 100000| 100|0.64332814| PASSED + diehard_runs| 0| 100000| 100|0.93746595| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.10e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_craps| 0| 200000| 100|0.00000000| FAILED + diehard_craps| 0| 200000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.12e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED + marsaglia_tsang_gcd| 0| 10000000| 100|0.00003253| WEAK +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.12e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_monobit| 1| 100000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.08e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_runs| 2| 100000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.11e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_serial| 1| 100000| 100|0.00000000| FAILED + sts_serial| 2| 100000| 100|0.00000000| FAILED + sts_serial| 3| 100000| 100|0.00000000| FAILED + sts_serial| 3| 100000| 100|0.00000000| FAILED + sts_serial| 4| 100000| 100|0.00000000| FAILED + sts_serial| 4| 100000| 100|0.00000000| FAILED + sts_serial| 5| 100000| 100|0.00000000| FAILED + sts_serial| 5| 100000| 100|0.00000000| FAILED + sts_serial| 6| 100000| 100|0.00000000| FAILED + sts_serial| 6| 100000| 100|0.00000000| FAILED + sts_serial| 7| 100000| 100|0.00000000| FAILED + sts_serial| 7| 100000| 100|0.00000000| FAILED + sts_serial| 8| 100000| 100|0.00000000| FAILED + sts_serial| 8| 100000| 100|0.00000000| FAILED + sts_serial| 9| 100000| 100|0.00000000| FAILED + sts_serial| 9| 100000| 100|0.00404845| WEAK + sts_serial| 10| 100000| 100|0.00000000| FAILED + sts_serial| 10| 100000| 100|0.00026880| WEAK + sts_serial| 11| 100000| 100|0.00000000| FAILED + sts_serial| 11| 100000| 100|0.00000112| WEAK + sts_serial| 12| 100000| 100|0.00000000| FAILED + sts_serial| 12| 100000| 100|0.00010025| WEAK + sts_serial| 13| 100000| 100|0.00000000| FAILED + sts_serial| 13| 100000| 100|0.00003788| WEAK + sts_serial| 14| 100000| 100|0.00000000| FAILED + sts_serial| 14| 100000| 100|0.00002685| WEAK + sts_serial| 15| 100000| 100|0.00000000| FAILED + sts_serial| 15| 100000| 100|0.00000295| WEAK + sts_serial| 16| 100000| 100|0.00000000| FAILED + sts_serial| 16| 100000| 100|0.00000071| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 6.10e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +rgb_minimum_distance| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 4.17e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_permutations| 5| 100000| 100|0.61242372| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.78e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_lagged_sum| 0| 1000000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.75e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_kstest_test| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.67e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_bytedistrib| 0| 51200000| 1|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.74e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_dct| 256| 50000| 1|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.78e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree| 32| 15000000| 1|0.35802483| PASSED + dab_filltree| 32| 15000000| 1|0.03798318| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.75e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree2| 0| 5000000| 1|0.00000000| FAILED + dab_filltree2| 1| 5000000| 1|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.0| 5.31e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_monobit2| 12| 65000000| 1|1.00000000| FAILED diff --git a/devtools/rngtest/results.1 b/devtools/rngtest/results.1 new file mode 100644 index 0000000000..1b93445751 --- /dev/null +++ b/devtools/rngtest/results.1 @@ -0,0 +1,260 @@ +Cpptraj type 1 +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.78e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_birthdays| 0| 100| 100|0.93575390| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.49e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_operm5| 0| 1000000| 100|0.45709550| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 4.96e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_32x32| 0| 40000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.82e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_6x8| 0| 100000| 100|0.97355690| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.83e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_str| 0| 256000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.85e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_byt| 0| 256000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_parking_lot| 0| 12000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.85e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_2dsphere| 2| 8000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_3dsphere| 3| 4000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_squeeze| 0| 100000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_runs| 0| 100000| 100|0.58542322| PASSED + diehard_runs| 0| 100000| 100|0.06198327| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.82e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_craps| 0| 200000| 100|0.00000000| FAILED + diehard_craps| 0| 200000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.83e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED + marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.77e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_monobit| 1| 100000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.80e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_runs| 2| 100000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.78e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_serial| 1| 100000| 100|0.00000000| FAILED + sts_serial| 2| 100000| 100|0.00000000| FAILED + sts_serial| 3| 100000| 100|0.00000000| FAILED + sts_serial| 3| 100000| 100|0.06655926| PASSED + sts_serial| 4| 100000| 100|0.00000000| FAILED + sts_serial| 4| 100000| 100|0.47071956| PASSED + sts_serial| 5| 100000| 100|0.00000000| FAILED + sts_serial| 5| 100000| 100|0.32512416| PASSED + sts_serial| 6| 100000| 100|0.00000000| FAILED + sts_serial| 6| 100000| 100|0.90811939| PASSED + sts_serial| 7| 100000| 100|0.00000000| FAILED + sts_serial| 7| 100000| 100|0.77041263| PASSED + sts_serial| 8| 100000| 100|0.00000000| FAILED + sts_serial| 8| 100000| 100|0.54120617| PASSED + sts_serial| 9| 100000| 100|0.00000000| FAILED + sts_serial| 9| 100000| 100|0.89253138| PASSED + sts_serial| 10| 100000| 100|0.00000000| FAILED + sts_serial| 10| 100000| 100|0.48285019| PASSED + sts_serial| 11| 100000| 100|0.00000000| FAILED + sts_serial| 11| 100000| 100|0.25756005| PASSED + sts_serial| 12| 100000| 100|0.00000000| FAILED + sts_serial| 12| 100000| 100|0.78336661| PASSED + sts_serial| 13| 100000| 100|0.00000000| FAILED + sts_serial| 13| 100000| 100|0.89063909| PASSED + sts_serial| 14| 100000| 100|0.00000000| FAILED + sts_serial| 14| 100000| 100|0.19063191| PASSED + sts_serial| 15| 100000| 100|0.00000000| FAILED + sts_serial| 15| 100000| 100|0.67736151| PASSED + sts_serial| 16| 100000| 100|0.00000000| FAILED + sts_serial| 16| 100000| 100|0.01145560| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.79e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +rgb_minimum_distance| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.81e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_permutations| 5| 100000| 100|0.60573210| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.77e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_lagged_sum| 0| 1000000| 100|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.78e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_kstest_test| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.75e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_bytedistrib| 0| 51200000| 1|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.80e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_dct| 256| 50000| 1|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.81e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree| 32| 15000000| 1|0.50234459| PASSED + dab_filltree| 32| 15000000| 1|0.59811692| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.75e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree2| 0| 5000000| 1|0.00000000| FAILED + dab_filltree2| 1| 5000000| 1|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.1| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_monobit2| 12| 65000000| 1|1.00000000| FAILED diff --git a/devtools/rngtest/results.2 b/devtools/rngtest/results.2 new file mode 100644 index 0000000000..d2ef63a931 --- /dev/null +++ b/devtools/rngtest/results.2 @@ -0,0 +1,260 @@ +Cpptraj type 2 +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.89e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_birthdays| 0| 100| 100|0.99126512| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.87e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_operm5| 0| 1000000| 100|0.33616728| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_32x32| 0| 40000| 100|0.61227793| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.87e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_6x8| 0| 100000| 100|0.24651799| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_str| 0| 256000| 100|0.36963644| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.87e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_byt| 0| 256000| 100|0.75489517| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.89e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_parking_lot| 0| 12000| 100|0.36905637| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_2dsphere| 2| 8000| 100|0.83220542| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_3dsphere| 3| 4000| 100|0.27072439| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_squeeze| 0| 100000| 100|0.84348080| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.83e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_runs| 0| 100000| 100|0.38180757| PASSED + diehard_runs| 0| 100000| 100|0.15389951| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_craps| 0| 200000| 100|0.02194641| PASSED + diehard_craps| 0| 200000| 100|0.56562021| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED + marsaglia_tsang_gcd| 0| 10000000| 100|0.00013155| WEAK +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.89e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_monobit| 1| 100000| 100|0.65973052| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_runs| 2| 100000| 100|0.20210136| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_serial| 1| 100000| 100|0.65973052| PASSED + sts_serial| 2| 100000| 100|0.39276488| PASSED + sts_serial| 3| 100000| 100|0.67449310| PASSED + sts_serial| 3| 100000| 100|0.97109028| PASSED + sts_serial| 4| 100000| 100|0.58078908| PASSED + sts_serial| 4| 100000| 100|0.93722471| PASSED + sts_serial| 5| 100000| 100|0.91169891| PASSED + sts_serial| 5| 100000| 100|0.79827105| PASSED + sts_serial| 6| 100000| 100|0.80955694| PASSED + sts_serial| 6| 100000| 100|0.98277956| PASSED + sts_serial| 7| 100000| 100|0.82852678| PASSED + sts_serial| 7| 100000| 100|0.60788851| PASSED + sts_serial| 8| 100000| 100|0.98814112| PASSED + sts_serial| 8| 100000| 100|0.54498283| PASSED + sts_serial| 9| 100000| 100|0.22907191| PASSED + sts_serial| 9| 100000| 100|0.24623161| PASSED + sts_serial| 10| 100000| 100|0.51963051| PASSED + sts_serial| 10| 100000| 100|0.86614025| PASSED + sts_serial| 11| 100000| 100|0.97525698| PASSED + sts_serial| 11| 100000| 100|0.74434975| PASSED + sts_serial| 12| 100000| 100|0.05242035| PASSED + sts_serial| 12| 100000| 100|0.69318146| PASSED + sts_serial| 13| 100000| 100|0.39024353| PASSED + sts_serial| 13| 100000| 100|0.48925517| PASSED + sts_serial| 14| 100000| 100|0.17514209| PASSED + sts_serial| 14| 100000| 100|0.83802927| PASSED + sts_serial| 15| 100000| 100|0.31763244| PASSED + sts_serial| 15| 100000| 100|0.47683051| PASSED + sts_serial| 16| 100000| 100|0.57382841| PASSED + sts_serial| 16| 100000| 100|0.69274018| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +rgb_minimum_distance| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_permutations| 5| 100000| 100|0.15940518| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.87e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_lagged_sum| 0| 1000000| 100|0.50417233| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.88e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_kstest_test| 0| 10000| 1000|0.73392878| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.87e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_bytedistrib| 0| 51200000| 1|0.96427985| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_dct| 256| 50000| 1|0.44065572| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.89e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree| 32| 15000000| 1|0.43781104| PASSED + dab_filltree| 32| 15000000| 1|0.11050760| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree2| 0| 5000000| 1|0.46915097| PASSED + dab_filltree2| 1| 5000000| 1|0.52770543| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.2| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_monobit2| 12| 65000000| 1|0.71206595| PASSED diff --git a/devtools/rngtest/results.3 b/devtools/rngtest/results.3 new file mode 100644 index 0000000000..3c6e1d211b --- /dev/null +++ b/devtools/rngtest/results.3 @@ -0,0 +1,260 @@ +Cpptraj type 3 +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.80e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_birthdays| 0| 100| 100|0.06587073| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_operm5| 0| 1000000| 100|0.40149705| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_32x32| 0| 40000| 100|0.67447876| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.88e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_6x8| 0| 100000| 100|0.78021623| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_str| 0| 256000| 100|0.96511462| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_byt| 0| 256000| 100|0.19733905| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_parking_lot| 0| 12000| 100|0.88713368| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.89e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_2dsphere| 2| 8000| 100|0.43965991| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_3dsphere| 3| 4000| 100|0.38987829| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_squeeze| 0| 100000| 100|0.94527996| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_runs| 0| 100000| 100|0.72897299| PASSED + diehard_runs| 0| 100000| 100|0.96384707| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_craps| 0| 200000| 100|0.13020621| PASSED + diehard_craps| 0| 200000| 100|0.69505758| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.88e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + marsaglia_tsang_gcd| 0| 10000000| 100|0.00026818| WEAK + marsaglia_tsang_gcd| 0| 10000000| 100|0.00000005| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_monobit| 1| 100000| 100|0.29583959| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_runs| 2| 100000| 100|0.16821220| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_serial| 1| 100000| 100|0.29583959| PASSED + sts_serial| 2| 100000| 100|0.57834967| PASSED + sts_serial| 3| 100000| 100|0.45214070| PASSED + sts_serial| 3| 100000| 100|0.98218470| PASSED + sts_serial| 4| 100000| 100|0.92219842| PASSED + sts_serial| 4| 100000| 100|0.81517658| PASSED + sts_serial| 5| 100000| 100|0.78993152| PASSED + sts_serial| 5| 100000| 100|0.05547932| PASSED + sts_serial| 6| 100000| 100|0.93608621| PASSED + sts_serial| 6| 100000| 100|0.72292385| PASSED + sts_serial| 7| 100000| 100|0.57056039| PASSED + sts_serial| 7| 100000| 100|0.88755734| PASSED + sts_serial| 8| 100000| 100|0.45359684| PASSED + sts_serial| 8| 100000| 100|0.93400894| PASSED + sts_serial| 9| 100000| 100|0.05555853| PASSED + sts_serial| 9| 100000| 100|0.09703108| PASSED + sts_serial| 10| 100000| 100|0.40833328| PASSED + sts_serial| 10| 100000| 100|0.24460928| PASSED + sts_serial| 11| 100000| 100|0.65244767| PASSED + sts_serial| 11| 100000| 100|0.81847484| PASSED + sts_serial| 12| 100000| 100|0.88915238| PASSED + sts_serial| 12| 100000| 100|0.70810740| PASSED + sts_serial| 13| 100000| 100|0.49583173| PASSED + sts_serial| 13| 100000| 100|0.38972586| PASSED + sts_serial| 14| 100000| 100|0.91190178| PASSED + sts_serial| 14| 100000| 100|0.34725131| PASSED + sts_serial| 15| 100000| 100|0.95828482| PASSED + sts_serial| 15| 100000| 100|0.99591800| WEAK + sts_serial| 16| 100000| 100|0.86929126| PASSED + sts_serial| 16| 100000| 100|0.58841564| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +rgb_minimum_distance| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_permutations| 5| 100000| 100|0.14574692| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_lagged_sum| 0| 1000000| 100|0.71901865| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.89e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_kstest_test| 0| 10000| 1000|0.50119988| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_bytedistrib| 0| 51200000| 1|0.38175742| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_dct| 256| 50000| 1|0.68964353| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree| 32| 15000000| 1|0.32794175| PASSED + dab_filltree| 32| 15000000| 1|0.85841654| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.79e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree2| 0| 5000000| 1|0.81423763| PASSED + dab_filltree2| 1| 5000000| 1|0.99111312| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.3| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_monobit2| 12| 65000000| 1|0.46243573| PASSED diff --git a/devtools/rngtest/results.4 b/devtools/rngtest/results.4 new file mode 100644 index 0000000000..fa0e6d9b3e --- /dev/null +++ b/devtools/rngtest/results.4 @@ -0,0 +1,260 @@ +Cpptraj type 4 +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.86e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_birthdays| 0| 100| 100|0.80125785| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_operm5| 0| 1000000| 100|0.71179546| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.73e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_32x32| 0| 40000| 100|0.64515320| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.74e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_rank_6x8| 0| 100000| 100|0.58463671| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.74e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_str| 0| 256000| 100|0.36497704| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.75e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +diehard_count_1s_byt| 0| 256000| 100|0.02738060| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.74e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_parking_lot| 0| 12000| 100|0.59549008| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.72e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_2dsphere| 2| 8000| 100|0.78503696| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.74e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_3dsphere| 3| 4000| 100|0.55202329| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.76e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_squeeze| 0| 100000| 100|0.15640654| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.76e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_runs| 0| 100000| 100|0.37324636| PASSED + diehard_runs| 0| 100000| 100|0.36681249| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.71e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + diehard_craps| 0| 200000| 100|0.87266568| PASSED + diehard_craps| 0| 200000| 100|0.95923227| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.75e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED + marsaglia_tsang_gcd| 0| 10000000| 100|0.00001885| WEAK +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.90e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_monobit| 1| 100000| 100|0.68739066| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_runs| 2| 100000| 100|0.67974661| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + sts_serial| 1| 100000| 100|0.68739066| PASSED + sts_serial| 2| 100000| 100|0.92434869| PASSED + sts_serial| 3| 100000| 100|0.47493716| PASSED + sts_serial| 3| 100000| 100|0.60511119| PASSED + sts_serial| 4| 100000| 100|0.79286917| PASSED + sts_serial| 4| 100000| 100|0.29818065| PASSED + sts_serial| 5| 100000| 100|0.98471104| PASSED + sts_serial| 5| 100000| 100|0.93338676| PASSED + sts_serial| 6| 100000| 100|0.99270095| PASSED + sts_serial| 6| 100000| 100|0.78981370| PASSED + sts_serial| 7| 100000| 100|0.86644689| PASSED + sts_serial| 7| 100000| 100|0.83590155| PASSED + sts_serial| 8| 100000| 100|0.49303039| PASSED + sts_serial| 8| 100000| 100|0.55651946| PASSED + sts_serial| 9| 100000| 100|0.22127976| PASSED + sts_serial| 9| 100000| 100|0.15099306| PASSED + sts_serial| 10| 100000| 100|0.74548037| PASSED + sts_serial| 10| 100000| 100|0.96939483| PASSED + sts_serial| 11| 100000| 100|0.45903245| PASSED + sts_serial| 11| 100000| 100|0.89646157| PASSED + sts_serial| 12| 100000| 100|0.25182668| PASSED + sts_serial| 12| 100000| 100|0.02367326| PASSED + sts_serial| 13| 100000| 100|0.06473590| PASSED + sts_serial| 13| 100000| 100|0.20846743| PASSED + sts_serial| 14| 100000| 100|0.04235961| PASSED + sts_serial| 14| 100000| 100|0.34291170| PASSED + sts_serial| 15| 100000| 100|0.70322411| PASSED + sts_serial| 15| 100000| 100|0.66104574| PASSED + sts_serial| 16| 100000| 100|0.28610367| PASSED + sts_serial| 16| 100000| 100|0.11066776| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# +rgb_minimum_distance| 0| 10000| 1000|0.00000000| FAILED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.93e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_permutations| 5| 100000| 100|0.47952957| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_lagged_sum| 0| 1000000| 100|0.39288477| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + rgb_kstest_test| 0| 10000| 1000|0.00278698| WEAK +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.84e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_bytedistrib| 0| 51200000| 1|0.39014627| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.91e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_dct| 256| 50000| 1|0.99633239| WEAK +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree| 32| 15000000| 1|0.71522845| PASSED + dab_filltree| 32| 15000000| 1|0.65385761| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_filltree2| 0| 5000000| 1|0.25955268| PASSED + dab_filltree2| 1| 5000000| 1|0.08741026| PASSED +#=============================================================================# +# dieharder version 3.31.1 Copyright 2003 Robert G. Brown # +#=============================================================================# + rng_name | filename |rands/second| + file_input| temp.4| 5.92e+06 | +#=============================================================================# + test_name |ntup| tsamples |psamples| p-value |Assessment +#=============================================================================# + dab_monobit2| 12| 65000000| 1|0.02868272| PASSED diff --git a/devtools/rngtest/rngtest.cpp b/devtools/rngtest/rngtest.cpp new file mode 100644 index 0000000000..4087f4ea5d --- /dev/null +++ b/devtools/rngtest/rngtest.cpp @@ -0,0 +1,158 @@ +#include +#include +#include +#include +#include +#include "../../src/Random.h" +#include "../../src/CpptrajStdio.h" + +static int DieHard_Tests(Random_Number const& RNG, int itype, int iseed, int icount) { + printf("#==================================================================\n"); + printf("# Cpptraj generator %s seed = %i \n", + Random_Number::CurrentDefaultRngStr(), iseed); + printf("#==================================================================\n"); + printf("type: d\n"); + printf("count: %i\n", icount); + printf("numbit: 32\n"); + +/* + double d_imax = (double)std::numeric_limits::max(); + for (int ic = 0; ic < icount; ic++) { + double rn = RNG.rn_gen() * d_imax; + int in = (int)rn; + printf("%10i\n", in); + }*/ + for (int ic = 0; ic < icount; ic++) { + printf("%10u\n", RNG.rn_num()); + } + + //printf("\n"); + return 0; +} + +static int Shuffle_Tests(Random_Number const& RNG, int itype, int iseed, int icount) { + printf("Shuffle tests (%s)\n", Random_Number::CurrentDefaultRngStr()); + printf("\tSeed : %i\n", iseed); + printf("\tCount : %i\n", icount); + printf("\tType : %i\n", itype); + + std::vector nums; + nums.reserve(icount); + for (int i = 0; i < icount; i++) + nums.push_back( i ); + + RNG.ShufflePoints( nums ); + printf("Shuffled points:"); + for (std::vector::const_iterator it = nums.begin(); it != nums.end(); ++it) + printf(" %i", *it); + printf("\n"); + return 0; +} + +static int Range_Tests(Random_Number const& RNG, int itype, int iseed, int icount, int imax) { + printf("Range tests (%s)\n", Random_Number::CurrentDefaultRngStr()); + printf("\tSeed : %i\n", iseed); + printf("\tCount : %i\n", icount); + printf("\tType : %i\n", itype); + + unsigned int min = 1; + unsigned int max = (unsigned int)imax; + unsigned int width = max - min + 1; + std::vector Counts( width, 0 ); + double avg = (double)icount / (double)width; + printf("\tAverage: %g\n", avg); + unsigned int oob = 0; + for (int i = 0; i < icount; i++) { + unsigned int rn = RNG.rn_num_interval(min, max); + long int idx = (long int)rn - (long int)min; + if (idx < 0 || idx >= (long int)Counts.size()) { + printf("Random # %u out of bounds (%li)\n", rn, idx); + oob++; + } else { + Counts[idx]++; + } + } + + printf("Final counts (%u out of bounds):\n", oob); + unsigned int zeroCount = 0; + unsigned int num = min; + double sum = 0; + double absSum = 0; + for (std::vector::const_iterator it = Counts.begin(); it != Counts.end(); ++it, ++num) + { + if (*it < 1) { + zeroCount++; + } else { + double delta = ((double)(*it)) - avg; + printf("\t%10u : %10u (%g)\n", num, *it, delta); + sum += delta; + if (delta < 0) delta = -delta; + absSum += delta; + } + } + if (zeroCount > 0) + printf("%u BINS HAVE NO POPULATION.\n", zeroCount); + double nonZeroCount = (double)(Counts.size() - zeroCount); + sum /= nonZeroCount; + absSum /= nonZeroCount; + printf("Avg delta = %g\n", sum); + printf("|Avg| delta = %g\n", absSum); + return 0; +} + +// Test the random number generators in cpptraj +int main(int argc, char** argv) { + SetWorldSilent(true); + Random_Number RNG; + + enum ModeType {DIEHARD = 0, SHUFFLE, RANGE }; + ModeType mode = DIEHARD; + int iseed = 0; + int icount = 0; + int itype = 0; + int imax = 10; + // Parse options + for (int iarg = 1; iarg < argc; iarg++) { + std::string arg( argv[iarg] ); + if (arg == "-S") { + iseed = atoi( argv[++iarg] ); + } else if (arg == "-t") { + icount = atoi( argv[++iarg] ); + } else if (arg == "-r") { + itype = atoi( argv[++iarg] ); + } else if (arg == "--max") { + imax = atoi( argv[++iarg] ); + } else if (arg == "--mode") { + std::string marg(argv[++iarg]); + if (marg == "diehard") + mode = DIEHARD; + else if (marg == "shuffle") + mode = SHUFFLE; + else if (marg == "range") + mode = RANGE; + else { + fprintf(stderr,"Error: Unrecognized mode: %s\n", marg.c_str()); + return 1; + } + } + } + + // Setup RNG + Random_Number::RngType rt = (Random_Number::RngType)itype; + Random_Number::SetDefaultRng(rt); + if (RNG.rn_set( iseed )) { + fprintf(stderr, "Error: Could not set up RNG.\n"); + return 1; + } + int err = 0; + + if (mode == DIEHARD) { + err = DieHard_Tests(RNG, itype, iseed, icount); + } else if (mode == SHUFFLE) { + err = Shuffle_Tests(RNG, itype, iseed, icount); + } else if (mode == RANGE) { + err = Range_Tests(RNG, itype, iseed, icount, imax); + } + + return err; +} diff --git a/devtools/testu01/CreateLinks.sh b/devtools/testu01/CreateLinks.sh new file mode 100755 index 0000000000..5dbdb366bf --- /dev/null +++ b/devtools/testu01/CreateLinks.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +for FILE in ../../src/RNG_MersenneTwister.cpp ../../src/RNG.cpp ../../src/CpptrajStdio.cpp ../../src/Random.cpp ../../src/RNG_Stdlib.cpp ../../src/RNG_Marsaglia.cpp ../../src/RNG_PCG32.cpp ../../src/RNG_Xoshiro128pp.cpp ../../src/xoshiro128plusplus.cpp ; do + ln -s $FILE . +done + diff --git a/devtools/testu01/Makefile b/devtools/testu01/Makefile new file mode 100644 index 0000000000..8f61dbbeee --- /dev/null +++ b/devtools/testu01/Makefile @@ -0,0 +1,29 @@ +include ../../config.h + +TESTU01_DIR=$(TESTU01_HOME) + +#AR= ar cqs + +#CXX=g++ +CXXFLAGS += -I../../src + +SOURCES=cpptraj_rng.cpp RNG_MersenneTwister.cpp RNG.cpp CpptrajStdio.cpp Random.cpp RNG_Stdlib.cpp RNG_Marsaglia.cpp RNG_PCG32.cpp RNG_Xoshiro128pp.cpp xoshiro128plusplus.cpp + +OBJECTS=$(SOURCES:.cpp=.o) + +all: a.out + +main.o: main.c + $(CC) -c -o main.o -Wall -O3 -I$(TESTU01_DIR)/include main.c + +#cpptraj.a: $(OBJECTS) +# -/bin/rm -f cpptraj.a +# $(AR) cpptraj.a $(OBJECTS) + +a.out: $(OBJECTS) main.o + $(CXX) -o a.out main.o $(OBJECTS) -L$(TESTU01_DIR)/lib -ltestu01 -lprobdist -lmylib -lm +# gcc -std=c99 -Wall -O3 -o a.out cpptraj.a -L$(TESTU01_DIR)/lib -ltestu01 -lprobdist -lmylib -lm + +clean: + -/bin/rm *.o a.out + diff --git a/devtools/testu01/README.md b/devtools/testu01/README.md new file mode 100644 index 0000000000..fe64fb7818 --- /dev/null +++ b/devtools/testu01/README.md @@ -0,0 +1,158 @@ +TestU01 Results +=============== + +These are results from the TestU01 test suite for the various RNGs in cpptraj. + +Marsaglia +========= +The Marsaglia generator didnt complete tests. + +C Stdlib +======== +``` + Version: TestU01 1.2.3 + Generator: CPPTRAJ + Number of statistics: 144 + Total CPU time: 00:32:50.33 + The following tests gave p-values outside [0.001, 0.9990]: + (eps means a value < 1.0e-300): + (eps1 means a value < 1.0e-15): + + Test p-value + ---------------------------------------------- + 1 SerialOver, t = 2 eps + 2 SerialOver, t = 4 eps + 3 CollisionOver, t = 2 eps + 5 CollisionOver, t = 4 eps + 7 CollisionOver, t = 8 eps + 9 CollisionOver, t = 20 eps + 11 BirthdaySpacings, t = 2 eps + 12 BirthdaySpacings, t = 3 eps + 13 BirthdaySpacings, t = 4 eps + 14 BirthdaySpacings, t = 7 eps + 18 ClosePairs NP, t = 2 2.0e-4 + 18 ClosePairs mNP, t = 2 4.6e-66 + 18 ClosePairs mNP1, t = 2 8.4e-70 + 18 ClosePairs NJumps, t = 2 eps + 19 ClosePairs NP, t = 3 1.2e-7 + 19 ClosePairs mNP, t = 3 1.2e-134 + 19 ClosePairs mNP1, t = 3 1.3e-142 + 19 ClosePairs NJumps, t = 3 eps + 20 ClosePairs NP, t = 7 7.7e-10 + 20 ClosePairs mNP, t = 7 1.8e-79 + 20 ClosePairs mNP1, t = 7 1.9e-240 + 20 ClosePairs NJumps, t = 7 eps + 23 SimpPoker, d = 16 eps + 25 SimpPoker, d = 64 eps + 26 SimpPoker, d = 64 eps + 27 CouponCollector, d = 4 eps + 29 CouponCollector, d = 16 eps + 30 CouponCollector, d = 16 6.2e-8 + 31 Gap, r = 0 eps + 32 Gap, r = 27 eps + 33 Gap, r = 0 eps + 34 Gap, r = 22 eps + 41 MaxOft, t = 5 eps + 41 MaxOft AD, t = 5 3.2e-157 + 42 MaxOft, t = 10 eps + 42 MaxOft AD, t = 10 1.8e-79 + 43 MaxOft, t = 20 eps + 43 MaxOft AD, t = 20 1 - eps1 + 44 MaxOft, t = 30 eps + 44 MaxOft AD, t = 30 1 - eps1 + 45 SampleProd, t = 10 1 - eps1 + 46 SampleProd, t = 30 1 - eps1 + 47 SampleMean eps + 48 SampleCorr 1 - eps1 + 49 AppearanceSpacings, r = 0 1 - eps1 + 51 WeightDistrib, r = 0 eps + 52 WeightDistrib, r = 8 eps + 53 WeightDistrib, r = 16 eps + 54 WeightDistrib, r = 24 eps + 55 SumCollector eps + 56 MatrixRank, 60 x 60 eps + 58 MatrixRank, 300 x 300 eps + 60 MatrixRank, 1200 x 1200 eps + 62 Savir2 eps + 65 RandomWalk1 H (L = 90) eps + 65 RandomWalk1 M (L = 90) eps + 65 RandomWalk1 J (L = 90) eps + 65 RandomWalk1 R (L = 90) eps + 65 RandomWalk1 C (L = 90) eps + 67 RandomWalk1 H (L = 1000) eps + 67 RandomWalk1 M (L = 1000) eps + 67 RandomWalk1 J (L = 1000) eps + 67 RandomWalk1 R (L = 1000) eps + 67 RandomWalk1 C (L = 1000) eps + 69 RandomWalk1 H (L = 10000) eps + 69 RandomWalk1 M (L = 10000) eps + 69 RandomWalk1 J (L = 10000) eps + 69 RandomWalk1 R (L = 10000) eps + 69 RandomWalk1 C (L = 10000) eps + 71 LinearComp, r = 0 1 - eps1 + 71 LinearComp, r = 0 1 - eps1 + 73 LempelZiv 1 - eps1 + 74 Fourier3, r = 0 eps + 76 LongestHeadRun, r = 0 eps + 76 LongestHeadRun, r = 0 1 - 8.1e-9 + 80 HammingWeight2, r = 0 eps + 82 HammingCorr, L = 30 eps + 83 HammingCorr, L = 300 eps + 84 HammingCorr, L = 1200 eps + 85 HammingIndep, L = 30 eps + 87 HammingIndep, L = 300 eps + 89 HammingIndep, L = 1200 eps + 91 Run of bits, r = 0 eps + 95 AutoCor, d = 30 1 - eps1 + ---------------------------------------------- + All other tests were passed +``` + +Mersenne Twister +================ +``` + Version: TestU01 1.2.3 + Generator: CPPTRAJ + Number of statistics: 144 + Total CPU time: 00:31:10.94 + The following tests gave p-values outside [0.001, 0.9990]: + (eps means a value < 1.0e-300): + (eps1 means a value < 1.0e-15): + + Test p-value + ---------------------------------------------- + 67 RandomWalk1 C (L = 1000) 6.7e-4 + 71 LinearComp, r = 0 1 - eps1 + 72 LinearComp, r = 29 1 - eps1 + ---------------------------------------------- + All other tests were passed +``` + +PCG32 +===== +``` + Version: TestU01 1.2.3 + Generator: CPPTRAJ + Number of statistics: 144 + Total CPU time: 00:27:29.09 + The following tests gave p-values outside [0.001, 0.9990]: + (eps means a value < 1.0e-300): + (eps1 means a value < 1.0e-15): + + Test p-value + ---------------------------------------------- + 44 MaxOft, t = 30 3.6e-4 + ---------------------------------------------- + All other tests were passed +``` + +Xoshiro128++ +============ +``` + Version: TestU01 1.2.3 + Generator: CPPTRAJ + Number of statistics: 144 + Total CPU time: 00:28:41.95 + + All tests were passed +``` diff --git a/devtools/testu01/Test.sh b/devtools/testu01/Test.sh new file mode 100755 index 0000000000..51c107883d --- /dev/null +++ b/devtools/testu01/Test.sh @@ -0,0 +1,12 @@ +#!/bin/bash + +make clean +make + +MODE=2 +TESTLIST='0 1 2 3 4' + +for TEST in $TESTLIST ; do + echo "Test rng $TEST" + ./a.out --mode $MODE -r $TEST > results.$MODE.$TEST +done diff --git a/devtools/testu01/cpptraj_rng.cpp b/devtools/testu01/cpptraj_rng.cpp new file mode 100644 index 0000000000..5fce3708eb --- /dev/null +++ b/devtools/testu01/cpptraj_rng.cpp @@ -0,0 +1,19 @@ +#include "cpptraj_rng.h" +#include "../../src/Random.h" + +void* get_cpptraj_rng(int itype, int iseed) { + Random_Number::SetDefaultRng((Random_Number::RngType)itype); + Random_Number* rng = new Random_Number(); + rng->rn_set( iseed ); + + return ( reinterpret_cast( rng ) ); +} + +void destroy_cpptraj_rng(void* rng) { + delete( reinterpret_cast( rng ) ); +} + +unsigned int num_cpptraj_rng(void* rngIn) { + Random_Number* rng = reinterpret_cast( rngIn ); + return rng->rn_num(); +} diff --git a/devtools/testu01/cpptraj_rng.h b/devtools/testu01/cpptraj_rng.h new file mode 100644 index 0000000000..86bc65f5a3 --- /dev/null +++ b/devtools/testu01/cpptraj_rng.h @@ -0,0 +1,17 @@ +#ifndef INC_CPPTRAJ_RNG_H +#define INC_CPPTRAJ_RNG_H + +#ifdef __cplusplus +extern "C" { +#endif +/// Call with type, seed +void* get_cpptraj_rng(int,int); + +void destroy_cpptraj_rng(void*); + +unsigned int num_cpptraj_rng(void* ); +#ifdef __cplusplus +} +#endif + +#endif diff --git a/devtools/testu01/main.c b/devtools/testu01/main.c new file mode 100644 index 0000000000..0f0bc1c5b9 --- /dev/null +++ b/devtools/testu01/main.c @@ -0,0 +1,74 @@ +// Test the RNGs in cpptraj with TestU01 +#include +#include +#include "TestU01.h" +#include "cpptraj_rng.h" + +static void* rngptr = 0; + +unsigned int test_rng(void) +{ + return num_cpptraj_rng( rngptr ); +} + +int argIs(const char* arg, const char* key) +{ + if (strcmp(arg, key) == 0) return 1; + return 0; +} + +int main(int argc, char** argv) { + int it = 0; + int iseed = 0; + int icount = 0; + int itype = 0; + int imax = 10; + int imode = 1; // 1 = small crush + // Parse options + for (int iarg = 1; iarg < argc; iarg++) { + if (argIs(argv[iarg], "-S")) { + iseed = atoi( argv[++iarg] ); + } else if (argIs(argv[iarg], "-t")) { + icount = atoi( argv[++iarg] ); + } else if (argIs(argv[iarg], "-r")) { + itype = atoi( argv[++iarg] ); + } else if (argIs(argv[iarg], "--max")) { + imax = atoi( argv[++iarg] ); + } else if (argIs(argv[iarg], "--mode")) { + imode = atoi( argv[++iarg] ); + }/* else if (arg == "--mode") { + std::string marg(argv[++iarg]); + if (marg == "diehard") + mode = DIEHARD; + else if (marg == "shuffle") + mode = SHUFFLE; + else if (marg == "range") + mode = RANGE; + else { + fprintf(stderr,"Error: Unrecognized mode: %s\n", marg.c_str()); + return 1; + } + }*/ + } + + // Setup RNG + rngptr = get_cpptraj_rng( itype, iseed ); + + if (imode == 0) { + for (it = 0; it < icount; it++) + printf("%10u\n", num_cpptraj_rng( rngptr )); + } else { + // Create TestU01 PRNG object for our generator + unif01_Gen* gen = unif01_CreateExternGenBits((char*)"CPPTRAJ", test_rng); + // Run the tests. + if (imode == 1) + bbattery_SmallCrush(gen); + else if (imode == 2) + bbattery_Crush(gen); + // Clean up. + unif01_DeleteExternGenBits(gen); + } + + destroy_cpptraj_rng( rngptr ); + return 0; +} diff --git a/doc/cpptraj.bib b/doc/cpptraj.bib index c8ce9dd259..41f55461c0 100644 --- a/doc/cpptraj.bib +++ b/doc/cpptraj.bib @@ -360,3 +360,57 @@ @Article{Sitkoff94 year = {1994} } +@techreport{Oneill2014, + author = "Melissa E. O'Neill", + title = "PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation", + institution = "Harvey Mudd College", + address = "Claremont, CA", + number = "HMC-CS-2014-0905", + year = "2014", + month = Sep, + xurl = "https://www.cs.hmc.edu/tr/hmc-cs-2014-0905.pdf", +} + +@article{Vigna2017, +author = {Vigna, Sebastiano}, +doi = {10.1016/j.cam.2016.11.006}, +issn = {03770427}, +journal = {Journal of Computational and Applied Mathematics}, +month = {may}, +pages = {175--181}, +title = {{Further scramblings of Marsaglia's xorshift generators}}, +url = {https://linkinghub.elsevier.com/retrieve/pii/S0377042716305301}, +volume = {315}, +year = {2017} +} + +@article{Roe2020, +author = {Roe, Daniel R. and Brooks, Bernard R.}, +doi = {10.1063/5.0013849}, +file = {:D$\backslash$:/DropboxFolders/Dropbox/Manuscripts/Equil{\_}2019/Revision/5.0013849.pdf:pdf}, +issn = {0021-9606}, +journal = {The Journal of Chemical Physics}, +month = {aug}, +number = {5}, +pages = {054123}, +publisher = {AIP Publishing, LLC}, +title = {{A protocol for preparing explicitly solvated systems for stable molecular dynamics simulations}}, +url = {https://doi.org/10.1063/5.0013849 http://aip.scitation.org/doi/10.1063/5.0013849}, +volume = {153}, +year = {2020} +} + +@article{Roe2021, +author = {Roe, Daniel R. and Brooks, Bernard R.}, +doi = {10.1016/j.jmgm.2021.107832}, +file = {:D$\backslash$:/Documents/Papers/Analysis/1-s2.0-S1093326321000012-main.pdf:pdf}, +issn = {10933263}, +journal = {Journal of Molecular Graphics and Modelling}, +pages = {107832}, +publisher = {Elsevier Ltd}, +title = {{Improving the Speed of Volumetric Density Map Generation via Cubic Spline Interpolation}}, +url = {https://doi.org/10.1016/j.jmgm.2021.107832}, +volume = {104}, +year = {2021} +} + diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 28a0a4e455..95a08bb503 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -7136,7 +7136,7 @@ The following general commands are available: \begin_layout Standard \align center \begin_inset Tabular - + @@ -7569,6 +7569,26 @@ Print data set to screen. \begin_inset Text +\begin_layout Plain Layout +random +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Change default random number generator, create random sets. +\end_layout + +\end_inset + + + + +\begin_inset Text + \begin_layout Plain Layout readdata \end_layout @@ -10234,6 +10254,105 @@ and the output would look like: 2 0.630 \end_layout +\begin_layout Subsection +random +\end_layout + +\begin_layout LyX-Code +random [setdefault {marsaglia|stdlib|mt|pcg32|xo128}] +\end_layout + +\begin_layout LyX-Code + [createset count <#> settype {int|float01} [seed <#>]] +\end_layout + +\begin_deeper +\begin_layout Description +setdefault If specified, change the default random number generator (RNG). +\end_layout + +\begin_deeper +\begin_layout Description +marsaglia Use the Marsaglia RNG that is used in the Amber MD programs sander/pme +md. +\end_layout + +\begin_layout Description +stdlib Use the C standard library RNG. +\end_layout + +\begin_layout Description +mt Use the C++11 implementation of the Mersenne twister (mt19937); only + available with C++11 support. +\end_layout + +\begin_layout Description +pcg32 Use the 32 bit version of the Permuted Congruential Generator. +\begin_inset CommandInset citation +LatexCommand citep +key "Oneill2014" +literal "true" + +\end_inset + + +\end_layout + +\begin_layout Description +xo128 Use the Xoshiro128++ RNG. +\begin_inset CommandInset citation +LatexCommand citep +key "Vigna2017" +literal "true" + +\end_inset + + +\end_layout + +\end_deeper +\begin_layout Description +createset If specified, create a 1D data set filled with random numbers + of the specified type. +\end_layout + +\begin_deeper +\begin_layout Description + Name of created set. +\end_layout + +\begin_layout Description +count +\begin_inset space ~ +\end_inset + +<#> The number of elements to put into the set. +\end_layout + +\begin_layout Description +settype +\begin_inset space ~ +\end_inset + +{int|float01} Type of numbers to use; integer or floating point between + 0 and 1. +\end_layout + +\begin_layout Description +seed +\begin_inset space ~ +\end_inset + +<#> Optional seed for the RNG. +\end_layout + +\end_deeper +\end_deeper +\begin_layout Standard +This command can be used to set the default random number generator used + in CPPTRAJ, and/or create a 1D data set filled with random values. +\end_layout + \begin_layout Subsection readdata \begin_inset CommandInset label @@ -36271,20 +36390,14 @@ buffer \begin_layout Standard The calculation is sped up by using cubic splines to interpolate the exponential function when calculating the Gaussians. - The details are in Roe & Brooks, -\begin_inset Quotes eld -\end_inset +\begin_inset CommandInset citation +LatexCommand citep +key "Roe2021" +literal "true" -Improving the Speed of Volumetric Density Map Generation via Cubic Spline - Interpolation -\begin_inset Quotes erd \end_inset -, J. - Mol. - Graph. - Model. - (2021). + \end_layout \begin_layout Subsection @@ -36684,7 +36797,7 @@ trajin \begin_layout Standard \align center \begin_inset Tabular - + @@ -37163,6 +37276,35 @@ Calculate Kullback-Leibler divergence between two data sets. \begin_inset Text +\begin_layout Plain Layout +evalplateau +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Evaluate whether the data in a 1D set has reached a single exponential plateau. +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout + +\end_layout + +\end_inset + + + + +\begin_inset Text + \begin_layout Plain Layout FFT \end_layout @@ -41244,6 +41386,233 @@ divergence ds1 ds2 Calculate Kullback-Leibler divergence between specified data sets. \end_layout +\begin_layout Subsection +evalplateau +\end_layout + +\begin_layout LyX-Code +evalplateau [name ] [tol ] [valacut ] +\end_layout + +\begin_layout LyX-Code + [initpct ] [finalpct ] +\end_layout + +\begin_layout LyX-Code + [chisqcut ] [slopecut ] [maxit ] +\end_layout + +\begin_layout LyX-Code + [out ] [resultsout ] [statsout ] +\end_layout + +\begin_layout LyX-Code + ... +\end_layout + +\begin_deeper +\begin_layout Description +name +\begin_inset space ~ +\end_inset + +] Name for output data sets. +\end_layout + +\begin_layout Description +tol +\begin_inset space ~ +\end_inset + + Curve fitting tolerance. + Default 0.00001. +\end_layout + +\begin_layout Description +valacut +\begin_inset space ~ +\end_inset + + ( +\begin_inset Quotes eld +\end_inset + +Value A cutoff +\begin_inset Quotes erd +\end_inset + +) Cutoff for last half average vs estimated long term value. + Default 0.01. +\end_layout + +\begin_layout Description +initpct +\begin_inset space ~ +\end_inset + + The initial percentage of data to use for the initial density guess. + Default 1%. +\end_layout + +\begin_layout Description +finalpct +\begin_inset space ~ +\end_inset + + The final percentatge of data to use for the final density guess. + Default 50%. +\end_layout + +\begin_layout Description +chisqcut +\begin_inset space ~ +\end_inset + + Curve fit chi-squared cutoff. + Defualt 0.5. +\end_layout + +\begin_layout Description +slopecut +\begin_inset space ~ +\end_inset + + Final slope of fitted curve cutoff. + Default 0.000001. +\end_layout + +\begin_layout Description +maxit +\begin_inset space ~ +\end_inset + + Maximum number of iterations to perform during curve fit. +\end_layout + +\begin_layout Description +out +\begin_inset space ~ +\end_inset + + File to write data and fitted curve to. +\end_layout + +\begin_layout Description +resultsout +\begin_inset space ~ +\end_inset + + File to write plateau results to. +\end_layout + +\begin_layout Description +statsout +\begin_inset space ~ +\end_inset + + File to write curve fitting stats to. +\end_layout + +\begin_layout Description + Data sets to evaluate plateau for. +\end_layout + +\begin_layout Standard +Data Sets Created +\end_layout + +\begin_layout Description +[A0] The A0 (initial density) values. +\end_layout + +\begin_layout Description +[A1] The A1 (rate constant) values. +\end_layout + +\begin_layout Description +[A2] The A2 (final density) values. +\end_layout + +\begin_layout Description +[OneA1] One over the A1 (rate constant) values. +\end_layout + +\begin_layout Description +[corr] Curve fit correlation. +\end_layout + +\begin_layout Description +[vala] Difference between last half average of data vs final density + (A2). +\end_layout + +\begin_layout Description +[chisq] Chi-squared of the curve fit. +\end_layout + +\begin_layout Description +[pltime] Plateau time (time at which all cutoffs satisfied). +\end_layout + +\begin_layout Description +[fslope] Final slope of fitted curve. +\end_layout + +\begin_layout Description +[name] Input set legend. +\end_layout + +\begin_layout Description +[result] Final result: yes, no, err (error). +\end_layout + +\end_deeper +\begin_layout Standard +Attempt to determine if data has +\begin_inset Quotes eld +\end_inset + +plateaued +\begin_inset Quotes erd +\end_inset + +, i.e. + stopped changing significantly by the end of the set. + Currently the defaults are set up to evaluate density data as part of the + system preparation protocol described by Roe & Brooks. +\begin_inset CommandInset citation +LatexCommand citep +key "Roe2020" +literal "true" + +\end_inset + +\end_layout + \begin_layout Subsection fft \end_layout diff --git a/src/Analysis_EvalPlateau.cpp b/src/Analysis_EvalPlateau.cpp index 5142c269ed..7563192760 100644 --- a/src/Analysis_EvalPlateau.cpp +++ b/src/Analysis_EvalPlateau.cpp @@ -16,7 +16,6 @@ Analysis_EvalPlateau::Analysis_EvalPlateau() : maxIt_(0), debug_(0) { - SetHidden(true); } /** Aspects for output data sets (data_[]). KEEP IN SYNC WITH OdataType */ diff --git a/src/Command.cpp b/src/Command.cpp index 83dcd13e34..5da4e312a3 100644 --- a/src/Command.cpp +++ b/src/Command.cpp @@ -36,6 +36,7 @@ #include "Exec_ViewRst.h" #include "Exec_Set.h" #include "Exec_Show.h" +#include "Exec_Random.h" // ----- SYSTEM ---------------------------------------------------------------- #include "Exec_System.h" // ----- COORDS ---------------------------------------------------------------- @@ -228,6 +229,7 @@ void Command::Init() { Command::AddCmd( new Exec_PrintData(), Cmd::EXE, 1, "printdata" ); Command::AddCmd( new Exec_QuietBlocks(), Cmd::EXE, 1, "quietblocks" ); Command::AddCmd( new Exec_Quit(), Cmd::EXE, 2, "exit", "quit" ); + Command::AddCmd( new Exec_Random(), Cmd::EXE, 2, "random", "rng" ); Command::AddCmd( new Exec_ReadData(), Cmd::EXE, 1, "readdata" ); Command::AddCmd( new Exec_ReadEnsembleData(),Cmd::EXE, 1, "readensembledata" ); Command::AddCmd( new Exec_ReadInput(), Cmd::EXE, 1, "readinput" ); @@ -389,7 +391,7 @@ void Command::Init() { Command::AddCmd( new Analysis_CurveFit(), Cmd::ANA, 1, "curvefit" ); Command::AddCmd( new Analysis_Matrix(), Cmd::ANA, 2, "diagmatrix", "matrix" ); Command::AddCmd( new Analysis_Divergence(), Cmd::ANA, 1, "divergence" ); - Command::AddCmd( new Analysis_EvalPlateau(), Cmd::ANA, 1, "evalplateau" ); // hidden + Command::AddCmd( new Analysis_EvalPlateau(), Cmd::ANA, 1, "evalplateau" ); Command::AddCmd( new Analysis_FFT(), Cmd::ANA, 1, "fft" ); Command::AddCmd( new Analysis_HausdorffDistance,Cmd::ANA,1,"hausdorff" ); Command::AddCmd( new Analysis_Hist(), Cmd::ANA, 2, "hist", "histogram" ); diff --git a/src/Cpptraj.cpp b/src/Cpptraj.cpp index 32dfcc494c..03f8844699 100644 --- a/src/Cpptraj.cpp +++ b/src/Cpptraj.cpp @@ -65,6 +65,7 @@ void Cpptraj::Usage() { " [-h | --help] [-V | --version] [--defines] [-debug <#>]\n" " [--interactive] [--log ] [-tl]\n" " [-ms ] [-mr ] [--mask ] [--resmask ]\n" + " [--rng %s]\n" "\t-p : * Load as a topology file.\n" "\t-i : * Read input from file .\n" "\t-y : * Read from trajectory file ; same as input 'trajin '.\n" @@ -87,8 +88,9 @@ void Cpptraj::Usage() { "\t-mr : Print selected residue numbers to STDOUT.\n" "\t--mask : Print detailed atom selection to STDOUT.\n" "\t--resmask : Print detailed residue selection to STDOUT.\n" + "\t--rng : Change default random number generator.\n" " * Denotes flag may be specified multiple times.\n" - "\n"); + "\n", CpptrajState::RngKeywords()); } void Cpptraj::Intro() const { @@ -232,6 +234,9 @@ std::string Cpptraj::Defines() { #ifdef FFTW_FFT defined_str.append(" -DFFTW_FFT"); #endif +#ifdef C11_SUPPORT + defined_str.append(" -DC11_SUPPORT"); +#endif #ifdef LIBPME defined_str.append(" -DLIBPME"); #endif @@ -479,6 +484,9 @@ Cpptraj::Mode Cpptraj::ProcessCmdLineArgs(int argc, char** argv) { // --resmask: Parse mask string, print selected residue details if (ProcessMask( topFiles, refFiles, cmdLineArgs[++iarg], true, true )) return ERROR; return QUIT; + } else if ( NotFinalArg(cmdLineArgs, "--rng", iarg) ) { + // --rng: Change default RNG + if (State_.ChangeDefaultRng( cmdLineArgs[++iarg] )) return ERROR; } else { // Check if this is a file. # ifdef MPI diff --git a/src/CpptrajState.cpp b/src/CpptrajState.cpp index 6901f1dff9..4ea50c111f 100644 --- a/src/CpptrajState.cpp +++ b/src/CpptrajState.cpp @@ -5,6 +5,7 @@ #include "DataSet_Topology.h" // AddTopology #include "FrameArray.h" // RunEnsemble #include "ProgressBar.h" +#include "Random.h" #ifdef MPI # include "Parallel.h" # include "DataSet_Coords_TRJ.h" @@ -30,6 +31,38 @@ CpptrajState::CpptrajState() : # endif {} +/** \return Keywords recognized by ChangeDefaultRng() */ +const char* CpptrajState::RngKeywords() { +# ifdef C11_SUPPORT + return "{marsaglia|stdlib|mt|pcg32|xo128}"; +# else + return "{marsaglia|stdlib|pcg32|xo128}"; +# endif +} + +/** Change default RNG. */ +int CpptrajState::ChangeDefaultRng(std::string const& setarg) const { + if (!setarg.empty()) { + if (setarg == "marsaglia") Random_Number::SetDefaultRng( Random_Number::MARSAGLIA ); + else if (setarg == "stdlib") Random_Number::SetDefaultRng( Random_Number::STDLIB ); + else if (setarg == "mt") { +# ifdef C11_SUPPORT + Random_Number::SetDefaultRng( Random_Number::MERSENNE_TWISTER ); +# else + mprinterr("Error: Mersenne twister RNG requires C++11 support.\n"); + return 1; +# endif + } else if (setarg == "pcg32") Random_Number::SetDefaultRng( Random_Number::PCG32 ); + else if (setarg == "xo128") Random_Number::SetDefaultRng( Random_Number::XOSHIRO128PP ); + else { + mprinterr("Error: Unrecognized RNG type: %s\n", setarg.c_str()); + return 1; + } + mprintf("\tDefault RNG set to '%s'\n", Random_Number::CurrentDefaultRngStr()); + } + return 0; +} + /** This version of SetTrajMode() may be called from AddToActionQueue() to set * NORMAL as the default mode without setting up a trajectory. May also be used * to clear the mode via UNDEFINED. It is illegal to call this for ENSEMBLE. diff --git a/src/CpptrajState.h b/src/CpptrajState.h index 327ff3e6ca..b1f676fc75 100644 --- a/src/CpptrajState.h +++ b/src/CpptrajState.h @@ -21,6 +21,10 @@ class CpptrajState { void SetNoProgress() { showProgress_ = false; } void SetQuietBlocks(bool b) { quietBlocks_ = b; } void SetActionSilence(bool b) { actionList_.SetSilent(b); } + /// \return Keywords recognized by ChangeDefaultRng + static const char* RngKeywords(); + /// Change the default RNG according to given keyword. + int ChangeDefaultRng(std::string const&) const; # ifdef MPI void SetForceParaEnsemble(bool b) { forceParallelEnsemble_ = b; } # endif diff --git a/src/DataSet.cpp b/src/DataSet.cpp index f054fcafd8..04a76d739d 100644 --- a/src/DataSet.cpp +++ b/src/DataSet.cpp @@ -31,7 +31,8 @@ const char* DataSet::Descriptions_[] = { "parameters", // PARAMETERS "tensor", // TENSOR "string variable", // STRINGVAR - "vector with scalar" // VECTOR_SCALAR + "vector with scalar", // VECTOR_SCALAR + "unsigned integer" // UNSIGNED_INTEGER }; // CONSTRUCTOR diff --git a/src/DataSet.h b/src/DataSet.h index 27b5aa57e2..4a411e930b 100644 --- a/src/DataSet.h +++ b/src/DataSet.h @@ -28,7 +28,7 @@ class DataSet { UNKNOWN_DATA=0, DOUBLE, FLOAT, INTEGER, STRING, MATRIX_DBL, MATRIX_FLT, COORDS, VECTOR, MODES, GRID_FLT, GRID_DBL, REMLOG, XYMESH, TRAJ, REF_FRAME, MAT3X3, TOPOLOGY, CMATRIX, CMATRIX_NOMEM, CMATRIX_DISK, PH, PH_EXPL, PH_IMPL, - PARAMETERS, TENSOR, STRINGVAR, VECTOR_SCALAR + PARAMETERS, TENSOR, STRINGVAR, VECTOR_SCALAR, UNSIGNED_INTEGER }; /// Group DataSet belongs to. enum DataGroup { diff --git a/src/DataSetList.cpp b/src/DataSetList.cpp index 323f4995eb..a8e613a384 100644 --- a/src/DataSetList.cpp +++ b/src/DataSetList.cpp @@ -34,6 +34,7 @@ #include "DataSet_Tensor.h" #include "DataSet_StringVar.h" #include "DataSet_Vector_Scalar.h" +#include "DataSet_unsignedInt.h" bool DataSetList::useDiskCache_ = false; @@ -59,6 +60,7 @@ DataSet* DataSetList::NewSet(DataSet::DataType typeIn) { ds = DataSet_integer_mem::Alloc(); # endif break; + case DataSet::UNSIGNED_INTEGER : ds = DataSet_unsignedInt::Alloc(); break; case DataSet::STRING : ds = DataSet_string::Alloc(); break; case DataSet::MATRIX_DBL : ds = DataSet_MatrixDbl::Alloc(); break; case DataSet::MATRIX_FLT : ds = DataSet_MatrixFlt::Alloc(); break; diff --git a/src/DataSet_unsignedInt.cpp b/src/DataSet_unsignedInt.cpp new file mode 100644 index 0000000000..4fb6cc0ac9 --- /dev/null +++ b/src/DataSet_unsignedInt.cpp @@ -0,0 +1,100 @@ +#include "DataSet_unsignedInt.h" +//#inc lude "CpptrajStdio.h" // DEBUG + +/** Reserve space in the Data and Frames arrays. */ +int DataSet_unsignedInt::Allocate( SizeArray const& sizeIn ) { + if (!sizeIn.empty()) + Data_.reserve( sizeIn[0] ); + return 0; +} + +/** Allocate space in Data_ array. */ +int DataSet_unsignedInt::MemAlloc( SizeArray const& sizeIn ) { + if (!sizeIn.empty()) + Data_.resize( sizeIn[0] ); + return 0; +} + +/** Copy a block of data of nelts elements from set dptrIn at position pos to startIdx. */ +void DataSet_unsignedInt::CopyBlock(size_t startIdx, DataSet const* dptrIn, size_t pos, size_t nelts) +{ + DataSet_unsignedInt const& setIn = static_cast( *dptrIn ); + const unsigned int* ptr = (&(setIn.Data_[0])+pos); + std::copy( ptr, ptr + nelts, &(Data_[0]) + startIdx ); +} + +/** Insert data vIn at frame. */ +void DataSet_unsignedInt::Add(size_t frame, const void* vIn) { + if (frame > Data_.size()) + Data_.resize( frame, 0 ); + // Always insert at the end + // NOTE: No check for duplicate frame values. + Data_.push_back( *((unsigned int*)vIn) ); +} + +/** Write data at frame to CharBuffer. If no data for frame write 0.0. + */ +void DataSet_unsignedInt::WriteBuffer(CpptrajFile &cbuffer, SizeArray const& pIn) const { + if (pIn[0] >= Data_.size()) + cbuffer.Printf(format_.fmt(), 0); + else + cbuffer.Printf(format_.fmt(), Data_[pIn[0]]); +} + +/** Append set of same kind. */ +int DataSet_unsignedInt::Append(DataSet* dsIn) { + if (dsIn->Empty()) return 0; + if (dsIn->Group() != SCALAR_1D) return 1; + if (dsIn->Type() == UNSIGNED_INTEGER) { + size_t oldsize = Size(); + std::vector const& dataIn = ((DataSet_unsignedInt*)dsIn)->Data_; + Data_.resize( oldsize + dataIn.size() ); + std::copy( dataIn.begin(), dataIn.end(), Data_.begin() + oldsize ); + } else { + DataSet_1D const& ds = static_cast( *dsIn ); + for (unsigned int i = 0; i != ds.Size(); i++) + Data_.push_back( (unsigned int)ds.Dval(i) ); + } + return 0; +} + +#ifdef MPI +// DataSet_unsignedInt::Sync() +int DataSet_unsignedInt::Sync(size_t total, std::vector const& rank_frames, + Parallel::Comm const& commIn) +{ + if (commIn.Size()==1) return 0; + if (commIn.Master()) { + // Resize for total number of frames. + Data_.resize( total ); + unsigned int* endptr = &(Data_[0]) + rank_frames[0]; + // Receive data from each rank + for (int rank = 1; rank < commIn.Size(); rank++) { + commIn.SendMaster( endptr, rank_frames[rank], rank, MPI_UNSIGNED ); + endptr += rank_frames[rank]; + } + } else // Send data to master //TODO adjust for repeated additions? + commIn.SendMaster( &(Data_[0]), Data_.size(), commIn.Rank(), MPI_UNSIGNED ); + return 0; +} + +/** Receive unsigned integer data from rank, append to after offset. */ +int DataSet_unsignedInt::Recv(size_t total, unsigned int offset, int nframes, + int fromRank, int tag, Parallel::Comm const& commIn) +{ + Data_.resize( total ); + unsigned int* d_beg = &(Data_[0]) + offset; + //rprintf("DEBUG: Receive %i frames fromRank %i tag %i\n", nframes, fromRank, tag); + if (commIn.Recv( d_beg, nframes, MPI_UNSIGNED, fromRank, tag )) return 1; + // TODO Do SetNeedsSync false here? + return 0; +} + +/** Send unsigned integer data to rank. */ +int DataSet_unsignedInt::Send(int toRank, int tag, Parallel::Comm const& commIn) +const +{ + //rprintf("DEBUG: Send %zu frames toRank %i tag %i\n", Data_.size(), toRank, tag); + return commIn.Send( (void*)&(Data_[0]), Data_.size(), MPI_UNSIGNED, toRank, tag ); +} +#endif diff --git a/src/DataSet_unsignedInt.h b/src/DataSet_unsignedInt.h new file mode 100644 index 0000000000..940d9c383a --- /dev/null +++ b/src/DataSet_unsignedInt.h @@ -0,0 +1,56 @@ +#ifndef INC_DATASET_UNSIGNEDINT_MEM_H +#define INC_DATASET_UNSIGNEDINT_MEM_H +#include +#include "DataSet_1D.h" +/// Hold an array of unsigned integer values in memory. +class DataSet_unsignedInt : public DataSet_1D { + public: + DataSet_unsignedInt() : DataSet_1D(UNSIGNED_INTEGER, TextFormat(TextFormat::UNSIGNED, 12)) {} + static DataSet* Alloc() { return (DataSet*)new DataSet_unsignedInt();} + // ----- DataSet_integer functions ----------- + void SetElement(size_t idx, unsigned int val) { Data_[idx] = val; } + unsigned int operator[](size_t idx) const { return Data_[idx]; } + void AddElement(unsigned int i) { Data_.push_back( i ); } + /// Make set size sizeIn, all values set to 0.0. + void Resize(size_t sizeIn) { Data_.resize(sizeIn, 0); } + /// Make set size sizeIn, all values set to val. + void Assign(size_t sizeIn, unsigned int val) { Data_.resize(sizeIn, val); } + inline void AddVal(size_t, unsigned int); + // ----- DataSet_1D functions ---------------- + double Xcrd(size_t idx) const { return Dim(0).Coord(idx); } + void SetY(size_t idx, double y) { SetElement(idx, (unsigned int)y); } + double Dval(size_t idx) const { return (double)Data_[idx]; } + const void* VoidPtr(size_t idx) const { return (void*)(&(Data_[0])+idx); } + // ----- DataSet functions ------------------- + size_t Size() const { return Data_.size(); } +# ifdef MPI + int Sync(size_t, std::vector const&, Parallel::Comm const&); + int Recv(size_t, unsigned int, int, int, int, Parallel::Comm const&); + int Send(int, int, Parallel::Comm const&) const; +# endif + void Info() const { return; } + int Allocate(SizeArray const&); + void Add( size_t, const void* ); + void WriteBuffer(CpptrajFile&, SizeArray const&) const; + int Append(DataSet*); + size_t MemUsageInBytes() const { return Data_.size() * sizeof(unsigned int); } + int MemAlloc(SizeArray const&); + void CopyBlock(size_t, const DataSet*, size_t, size_t); + // ------------------------------------------- + //typedef std::vector::iterator iterator; + //iterator begin() { return Data_.begin(); } + //iterator end() { return Data_.end(); } + //int* Ptr() { return &(Data_[0]); } + private: + std::vector Data_; +}; +// ----- INLINE FUNCTIONS ------------------------------------------------------ +void DataSet_unsignedInt::AddVal(size_t frame, unsigned int ival) { + if (frame < Data_.size()) + Data_[frame] += ival; + else { + if (frame > Data_.size()) Data_.resize( frame, 0 ); + Data_.push_back( ival ); + } +} +#endif diff --git a/src/Exec_Random.cpp b/src/Exec_Random.cpp new file mode 100644 index 0000000000..e322633396 --- /dev/null +++ b/src/Exec_Random.cpp @@ -0,0 +1,62 @@ +#include "Exec_Random.h" +#include "CpptrajStdio.h" +#include "Random.h" + +// Exec_Random::Help() +void Exec_Random::Help() const +{ + mprintf("\t[setdefault %s]\n", CpptrajState::RngKeywords()); + mprintf("\t[createset count <#> settype {int|float01} [seed <#>]]\n"); + mprintf(" If 'setdefault' specified, change the default random number generator used.\n" + " If 'createset' specified, create a 1D data set filled with random numbers\n" + " of the specified type: 'int' creates a set of integer numbers, 'float01'\n" + " creates a set of floating point numbers between 0 and 1.\n"); +} + +// Exec_Random::Execute() +Exec::RetType Exec_Random::Execute(CpptrajState& State, ArgList& argIn) +{ + std::string setarg = argIn.GetStringKey("setdefault"); + if (!setarg.empty()) { + if (State.ChangeDefaultRng( setarg )) return CpptrajState::ERR; + } + + std::string dsname = argIn.GetStringKey("createset"); + if (!dsname.empty()) { + int iseed = argIn.getKeyInt("seed", -1); + int count = argIn.getKeyInt("count", -1); + if (count < 1) { + mprinterr("Error: Must specify 'count' > 0 for 'createset'\n"); + return CpptrajState::ERR; + } + DataFile* outfile = State.DFL().AddDataFile( argIn.GetStringKey("out"), argIn ); + + Random_Number rng; + if (rng.rn_set( iseed )) return CpptrajState::ERR; + std::string typestr = argIn.GetStringKey("settype"); + if (typestr == "int") { + // Create integer set + DataSet* ds = State.DSL().AddSet( DataSet::UNSIGNED_INTEGER, MetaData(dsname) ); + if (ds == 0) return CpptrajState::ERR; + if (outfile != 0) outfile->AddDataSet( ds ); + for (int idx = 0; idx != count; idx++) { + unsigned int rn = rng.rn_num(); + ds->Add(idx, &rn); + } + } else if (typestr == "float01") { + // Create floating point set between 0 and 1 + DataSet* ds = State.DSL().AddSet( DataSet::DOUBLE, MetaData(dsname) ); + if (ds == 0) return CpptrajState::ERR; + if (outfile != 0) outfile->AddDataSet( ds ); + for (int idx = 0; idx != count; idx++) { + double rn = rng.rn_gen(); + ds->Add(idx, &rn); + } + } else { + mprinterr("Error: Unrecognized 'settype' for 'createset': %s\n", typestr.c_str()); + return CpptrajState::ERR; + } + } + + return CpptrajState::OK; +} diff --git a/src/Exec_Random.h b/src/Exec_Random.h new file mode 100644 index 0000000000..16863cfa1c --- /dev/null +++ b/src/Exec_Random.h @@ -0,0 +1,12 @@ +#ifndef INC_EXEC_RANDOM_H +#define INC_EXEC_RANDOM_H +#include "Exec.h" +/// Set RNG type; create sets with random data. +class Exec_Random : public Exec { + public: + Exec_Random() : Exec(GENERAL) {} + void Help() const; + DispatchObject* Alloc() const { return (DispatchObject*)new Exec_Random(); } + RetType Execute(CpptrajState&, ArgList&); +}; +#endif diff --git a/src/ExternalFxn.cpp b/src/ExternalFxn.cpp new file mode 100644 index 0000000000..f55948e8cb --- /dev/null +++ b/src/ExternalFxn.cpp @@ -0,0 +1,14 @@ +#include "Random.h" + +/** Allow external programs to change the default random number generator. + * Wrapping it here avoids e.g. pytraj having to use namespace Cpptraj + * which will fail since I made both a Cpptraj namespace and class. + * Its primary function right now is to allow pytraj tests to set + * the default RNG back to Marsaglia. + * It also avoids having to import the Random_Number::RngType enum. + * TODO this can be deprecated when the Cpptraj class is renamed. + */ +void EXT_SetDefaultRng(int rtypeIn) { + Random_Number::RngType rt = (Random_Number::RngType)rtypeIn; + Random_Number::SetDefaultRng( rt ); +} diff --git a/src/ExternalFxn.h b/src/ExternalFxn.h new file mode 100644 index 0000000000..a582b84524 --- /dev/null +++ b/src/ExternalFxn.h @@ -0,0 +1,8 @@ +#ifndef INC_EXTERNALFXN_H +#define INC_EXTERNALFXN_H +/*! \file ExternalFxn.h + \brief This file contains some functions that external programs (like pytraj) need. + */ +/// Calls Random_Number::SetDefaultRng() +void EXT_SetDefaultRng(int); +#endif diff --git a/src/RNG.cpp b/src/RNG.cpp new file mode 100644 index 0000000000..2ec4b041fe --- /dev/null +++ b/src/RNG.cpp @@ -0,0 +1,23 @@ +#include // clock +#include "RNG.h" +#include "CpptrajStdio.h" + +/** CONSTRUCTOR */ +Cpptraj::RNG::RNG() : iseed_(-1) {} + +/** Set up the seed, then perform any specific setup required by the RNG. */ +int Cpptraj::RNG::Set_Seed(int iseedIn) { + // If iseed <= 0 set the seed based on wallclock time + if (iseedIn <= 0 ) { + clock_t wallclock = clock(); + iseed_ = (int)wallclock; + mprintf("Random_Number: seed is <= 0, using wallclock time as seed (%i)\n", iseed_); + } else + iseed_ = iseedIn; + // Set up specific to RNG + if (SetupRng()) { + mprinterr("Error: RNG setup failed.\n"); + return 1; + } + return 0; +} diff --git a/src/RNG.h b/src/RNG.h new file mode 100644 index 0000000000..f6888fbae2 --- /dev/null +++ b/src/RNG.h @@ -0,0 +1,35 @@ +#ifndef INC_RNG_H +#define INC_RNG_H +namespace Cpptraj { +/// Interface to random number generators +class RNG { + public: + /// CONSTRUCTOR + RNG(); + /// DESTRUCTOR - virtual since inherited + virtual ~RNG() {} + /// Generate a random number between 0.0 and 1.0 + virtual double Generate() = 0; + /// Generate a random integer + virtual unsigned int Number() = 0; + + /// Initialize the random number generator with the given seed + int Set_Seed(int); + /// Initialize RN generator with 71277 (Amber default) + //void Set_Seed() { Set_Seed(71277); } + + /// Reorder elements in input array using Fisher-Yates shuffle + //void Fisher_Yates_Shuffle(std::vector&); + + /// \return true if RN generator has been set up. + bool IsSet() const { return (iseed_ != -1); } + /// \return Value of RN generator seed + int Seed() const { return iseed_; } + protected: + /// Any internal setup that needs to be done by the RNG + virtual int SetupRng() = 0; + private: + int iseed_; ///< The random number generator seed. -1 means use time +}; +} +#endif diff --git a/src/RNG_Marsaglia.cpp b/src/RNG_Marsaglia.cpp new file mode 100644 index 0000000000..31364ff9f6 --- /dev/null +++ b/src/RNG_Marsaglia.cpp @@ -0,0 +1,133 @@ +#include // div +#include +#include "RNG_Marsaglia.h" +#include "CpptrajStdio.h" + +// CONSTRUCTOR +/// Only sets setup state to false. Actual setup is rn_set +Cpptraj::RNG_Marsaglia::RNG_Marsaglia() { + RN_generator.iseed = -1; +} + +/** Max value for this RNG, 2^24 bits */ +const double Cpptraj::RNG_Marsaglia::rng_max_ = 16777216.0; + +/** Initialization routine for Marsaglias random number generator + * as implemented in Amber 3.0 Rev A by George Seibel. See doc in amrand. + * + * Testing: Call amrset with iseed = 54185253. This should result + * in is1 = 1802 and is2 = 9373. Call amrand 20000 times, then six + * more times, printing the six random numbers * 2**24 (ie, 4096*4096) + * They should be: (6f12.1) + * 6533892.0 14220222.0 7275067.0 6172232.0 8354498.0 10633180.0 + */ +int Cpptraj::RNG_Marsaglia::SetupRng() { + RN_generator.iseed = Seed(); + // Two internal seeds used in Marsaglia algorithm: + int is1, is2; + const int is1max = 31328; + const int is2max = 30081; + + div_t divresult; + + // Construct two internal seeds from single unbound seed. + // is1 and is2 are quotient and remainder of iseed/IS2MAX. + // Add one to keep zero and one results from both mapping to one. + divresult = div( RN_generator.iseed, is2max ); + divresult.quot++; + divresult.rem++; + // Ensure range: 1 <= is1 <= is1max + if ( divresult.quot > 1 ) + is1 = divresult.quot; + else + is1 = 1; + if ( is1max < is1 ) + is1 = is1max; + // Ensure range: 1 <= is2 <= is2max + if ( divresult.rem > 1 ) + is2 = divresult.rem; + else + is2 = 1; + if ( is2max < is2 ) + is2 = is2max; + + int i = ((is1/177) % 177) + 2; + int j = (is1 % 177) + 2; + int k = ((is2/169) % 178) + 1; + int l = (is2 % 169); + + for (int ii = 0; ii < 97; ii++) { + double s = 0.0; + double t = 0.5; + for (int jj = 0; jj < 24; jj++) { + //m = mod(mod(i*j, 179)*k, 179) + int mtemp = (i*j) % 179; + int m = (mtemp*k) % 179; + i = j; + j = k; + k = m; + l = (((53*l) + 1) % 169); + //if (mod(l*m, 64) .ge. 32) s = s + t + if ( ((l*m) % 64) >= 32 ) s += t; + //t = 0.5d0 * t + t *= 0.5; + } + RN_generator.u[ii] = s; + } + + RN_generator.c = 362436.0 / rng_max_; + RN_generator.cd = 7654321.0 / rng_max_; + RN_generator.cm = 16777213.0 / rng_max_; + + RN_generator.i97 = 96; + RN_generator.j97 = 32; + return 0; +} + +/** \author Portable Random number generator by George Marsaglia + * Original author: Amber 3.0 Rev A implementation by George Seibel + * + * This random number generator originally appeared in 'Toward a Universal + * Random Number Generator' by George Marsaglia and Arif Zaman. Florida + * State University Report: FSU-SCRI-87-50 (1987) + * + * It was later modified by F. James and published in 'A Review of Pseudo- + * random Number Generators' + * + * This is claimed to be the best known random number generator available. + * It passes ALL of the tests for random number generators and has a + * period of 2^144, is completely portable (gives bit identical results on + * all machines with at least 24-bit mantissas in the floating point + * representation). + * + * The algorithm is a combination of a Fibonacci sequence (with lags of 97 + * and 33, and operation "subtraction plus one, modulo one") and an + * "arithmetic sequence" (using subtraction). + * + * \return A random number between 0.0 and 1.0 + * \return -1.0 if the random number generator is not initialized. + */ +double Cpptraj::RNG_Marsaglia::Generate() { + if (!IsSet()) { + mprinterr("Error: random number generator not initialized."); + return -1.0; + } + + double uni = RN_generator.u[RN_generator.i97] - RN_generator.u[RN_generator.j97]; + if (uni < 0.0) uni += 1.0; + RN_generator.u[RN_generator.i97] = uni; + RN_generator.i97--; + if (RN_generator.i97 == -1) RN_generator.i97 = 96; + RN_generator.j97--; + if (RN_generator.j97 == -1) RN_generator.j97 = 96; + RN_generator.c -= RN_generator.cd; + if (RN_generator.c < 0.0) RN_generator.c += RN_generator.cm; + uni -= RN_generator.c; + if (uni < 0.0) uni += 1.0; + return uni; +} + +unsigned int Cpptraj::RNG_Marsaglia::Number() { + double rn = Generate(); + return (unsigned int)(rn * rng_max_); +} diff --git a/src/RNG_Marsaglia.h b/src/RNG_Marsaglia.h new file mode 100644 index 0000000000..7a70bf63e4 --- /dev/null +++ b/src/RNG_Marsaglia.h @@ -0,0 +1,57 @@ +#ifndef INC_RNG_MARSAGLIA_H +#define INC_RNG_MARSAGLIA_H +#include "RNG.h" +namespace Cpptraj { +/// Marsaglia's random number generator as implemented in Amber 3.0 Rev A +/** + * Adapted from fortran code in $AMBERHOME/src/pmemd/src/random.F90 + * + * This random number generator originally appeared in "Toward a Universal + * Random Number Generator" by George Marsaglia and Arif Zaman. Florida + * State University Report: FSU-SCRI-87-50 (1987) + * + * It was later modified by F. James and published in "A Review of Pseudo- + * random Number Generators" + * + * This is claimed to be the best known random number generator available. + * It passes ALL of the tests for random number generators and has a + * period of 2^144, is completely portable (gives bit identical results on + * all machines with at least 24-bit mantissas in the floating point + * representation). + * + * The algorithm is a combination of a Fibonacci sequence (with lags of 97 + * and 33, and operation "subtraction plus one, modulo one") and an + * "arithmetic sequence" (using subtraction). + */ +class RNG_Marsaglia : public RNG { + public: + RNG_Marsaglia(); + /// Generate a random number between 0.0 and 1.0 + double Generate(); + unsigned int Number(); + private: + /// The max for this RNG, 2^24 bits + static const double rng_max_; + /// Setup RNG + int SetupRng(); + /// Variables necessary for Marsaglia random number stream. + /** This is placed in a struct in case the state of the random number + * generator ever needs to be stored. + */ + struct random_state { + // Real variables in Marsaglia algorithm + double u[97]; + double c; + double cd; + double cm; + // Pointers into u() in Marsaglia algorithm + int i97; + int j97; + // Random seed; if -1 random generator has not been set + int iseed; + }; + /// Hold the state of the random number generator + random_state RN_generator; +}; +} +#endif diff --git a/src/RNG_MersenneTwister.cpp b/src/RNG_MersenneTwister.cpp new file mode 100644 index 0000000000..3e483b90b0 --- /dev/null +++ b/src/RNG_MersenneTwister.cpp @@ -0,0 +1,32 @@ +#include "RNG_MersenneTwister.h" +#include "CpptrajStdio.h" + +Cpptraj::RNG_MersenneTwister::RNG_MersenneTwister() +{} + +double Cpptraj::RNG_MersenneTwister::Generate() { +# ifdef C11_SUPPORT + unsigned urn = gen_(); + return (double)urn / (double)gen_.max(); +# else + return 0; +# endif +} + +unsigned int Cpptraj::RNG_MersenneTwister::Number() { +# ifdef C11_SUPPORT + return gen_(); +# else + return 0; +# endif +} + +int Cpptraj::RNG_MersenneTwister::SetupRng() { +# ifdef C11_SUPPORT + gen_.seed( Seed() ); + return 0; +# else + mprinterr("Error: Cpptraj was compiled without C++11 support. Cannot use Mersenne twister.\n"); + return 1; +# endif +} diff --git a/src/RNG_MersenneTwister.h b/src/RNG_MersenneTwister.h new file mode 100644 index 0000000000..48831b4d49 --- /dev/null +++ b/src/RNG_MersenneTwister.h @@ -0,0 +1,21 @@ +#ifndef INC_RNG_MERSENNETWISTER_H +#define INC_RNG_MERSENNETWISTER_H +#include "RNG.h" +#ifdef C11_SUPPORT +#include +#endif +namespace Cpptraj { +class RNG_MersenneTwister : public RNG { + public: + RNG_MersenneTwister(); + + double Generate(); + unsigned int Number(); + private: + int SetupRng(); +# ifdef C11_SUPPORT + std::mt19937 gen_; ///< The Mersenne twister engine +# endif +}; +} +#endif diff --git a/src/RNG_PCG32.cpp b/src/RNG_PCG32.cpp new file mode 100644 index 0000000000..08166e0b1d --- /dev/null +++ b/src/RNG_PCG32.cpp @@ -0,0 +1,29 @@ +#include "RNG_PCG32.h" +#include "CpptrajStdio.h" + +int Cpptraj::RNG_PCG32::SetupRng() { +# ifdef _MSC_VER + mprinterr("Error: PCG32 generator does not yet work under Visual Studio\n"); + return 1; +# else + rng_ = pcg32( Seed() ); + return 0; +# endif +} + +double Cpptraj::RNG_PCG32::Generate() { +# ifdef _MSC_VER + return 0; +# else + unsigned int rn = rng_(); + return (double)rn / (double)rng_.max(); +# endif +} + +unsigned int Cpptraj::RNG_PCG32::Number() { +# ifdef _MSC_VER + return 0; +# else + return rng_(); +# endif +} diff --git a/src/RNG_PCG32.h b/src/RNG_PCG32.h new file mode 100644 index 0000000000..8c8e7ba952 --- /dev/null +++ b/src/RNG_PCG32.h @@ -0,0 +1,22 @@ +#ifndef INC_RNG_PCG32_H +#define INC_RNG_PCG32_H +#include "RNG.h" +#ifndef _MSC_VER +#include "pcg_random.hpp" +#endif +namespace Cpptraj { + +class RNG_PCG32 : public RNG { + public: + RNG_PCG32() {} + + double Generate(); + unsigned int Number(); + private: + int SetupRng(); +# ifndef _MSC_VER + pcg32 rng_; +# endif +}; +} +#endif diff --git a/src/RNG_Stdlib.cpp b/src/RNG_Stdlib.cpp new file mode 100644 index 0000000000..4e1ab3fc41 --- /dev/null +++ b/src/RNG_Stdlib.cpp @@ -0,0 +1,21 @@ +#include "RNG_Stdlib.h" +#include + +Cpptraj::RNG_Stdlib::RNG_Stdlib() {} + +int Cpptraj::RNG_Stdlib::SetupRng() { + srand( Seed() ); + return 0; +} + +double Cpptraj::RNG_Stdlib::Generate() { + const double drand_max = (double)RAND_MAX; + int rn = rand(); + double drn = (double)rn; + + return (drn / drand_max); +} + +unsigned int Cpptraj::RNG_Stdlib::Number() { + return rand(); +} diff --git a/src/RNG_Stdlib.h b/src/RNG_Stdlib.h new file mode 100644 index 0000000000..03a614b455 --- /dev/null +++ b/src/RNG_Stdlib.h @@ -0,0 +1,16 @@ +#ifndef INC_RNG_STDLIB_H +#define INC_RNG_STDLIB_H +#include "RNG.h" +namespace Cpptraj { +/// RNG from stdlib +class RNG_Stdlib : public RNG { + public: + RNG_Stdlib(); + + double Generate(); + unsigned int Number(); + private: + int SetupRng(); +}; +} +#endif diff --git a/src/RNG_Xoshiro128pp.cpp b/src/RNG_Xoshiro128pp.cpp new file mode 100644 index 0000000000..717dbcfa1f --- /dev/null +++ b/src/RNG_Xoshiro128pp.cpp @@ -0,0 +1,17 @@ +#include "xoshiro128plusplus.h" +#include "RNG_Xoshiro128pp.h" +#include + +int Cpptraj::RNG_Xoshiro128pp::SetupRng() { + x128pp_seed( Seed() ); + return 0; +} + +unsigned int Cpptraj::RNG_Xoshiro128pp::Number() { + return x128pp_next(); +} + +double Cpptraj::RNG_Xoshiro128pp::Generate() { + unsigned int rn = Number(); + return (double)rn / (double)std::numeric_limits::max(); +} diff --git a/src/RNG_Xoshiro128pp.h b/src/RNG_Xoshiro128pp.h new file mode 100644 index 0000000000..88a5776c7e --- /dev/null +++ b/src/RNG_Xoshiro128pp.h @@ -0,0 +1,16 @@ +#ifndef INC_RNG_XOSHIRO128PP_H +#define INC_RNG_XOSHIRO128PP_H +#include "RNG.h" +namespace Cpptraj { + +class RNG_Xoshiro128pp : public RNG { + public: + RNG_Xoshiro128pp() {} + + double Generate(); + unsigned int Number(); + private: + int SetupRng(); +}; +} +#endif diff --git a/src/Random.cpp b/src/Random.cpp index c4039d646a..7c5b474fde 100644 --- a/src/Random.cpp +++ b/src/Random.cpp @@ -1,154 +1,138 @@ -// Random_Number -#include // clock -#include // div -#include // sqrt, log -#include "Random.h" +#include //sqrt, log #include "CpptrajStdio.h" +#include "Random.h" +#include "RNG.h" +#include "RNG_Marsaglia.h" +#include "RNG_Stdlib.h" +#include "RNG_MersenneTwister.h" +#include "RNG_PCG32.h" +#include "RNG_Xoshiro128pp.h" + +/** Starting default type. */ +Random_Number::RngType Random_Number::defaultType_ = Random_Number::MARSAGLIA; +// TODO test these, then enable. Should be better RNGs in the long run +//#ifdef C11_SUPPORT +// Random_Number::MERSENNE_TWISTER; +//#else +// Random_Number::STDLIB; +//#endif + +/** Starting default seed. */ +//int Random_Number::defaultSeed_ = 71277; // AMBER default seed + +/** Set the default RNG type. */ +void Random_Number::SetDefaultRng(RngType r) { + defaultType_ = r; +} + +/** Set the default seed. */ +//void Random_Number::SetDefaultSeed(int i) { +// defaultSeed_ = i; +//} -// CONSTRUCTOR -/// Only sets setup state to false. Actual setup is rn_set -Random_Number::Random_Number() { - RN_generator.iseed = -1; +/** CONSTRUCTOR */ +Random_Number::Random_Number() : + rng_(0) +{} + +/** DESTRUCTOR */ +Random_Number::~Random_Number() { + if (rng_ != 0) delete rng_; } -// Random_Number::rn_set() -/** Initialization routine for Marsaglias random number generator - * as implemented in Amber 3.0 Rev A by George Seibel. See doc in amrand. - * - * Testing: Call amrset with iseed = 54185253. This should result - * in is1 = 1802 and is2 = 9373. Call amrand 20000 times, then six - * more times, printing the six random numbers * 2**24 (ie, 4096*4096) - * They should be: (6f12.1) - * 6533892.0 14220222.0 7275067.0 6172232.0 8354498.0 10633180.0 - * \param iseed Seed for the random number generator. - */ -void Random_Number::rn_set(int iseedIn) { - // Two internal seeds used in Marsaglia algorithm: - int is1, is2; - const int is1max = 31328; - const int is2max = 30081; - - div_t divresult; - - - // If iseed <= 0 set the seed based on wallclock time - if (iseedIn <= 0 ) { - clock_t wallclock = clock(); - RN_generator.iseed = (int) wallclock; - mprintf("Random_Number: seed is <= 0, using wallclock time as seed (%i)\n",RN_generator.iseed); - } else - RN_generator.iseed = iseedIn; - - // Construct two internal seeds from single unbound seed. - // is1 and is2 are quotient and remainder of iseed/IS2MAX. - // Add one to keep zero and one results from both mapping to one. - divresult = div( RN_generator.iseed, is2max ); - divresult.quot++; - divresult.rem++; - // Ensure range: 1 <= is1 <= is1max - if ( divresult.quot > 1 ) - is1 = divresult.quot; - else - is1 = 1; - if ( is1max < is1 ) - is1 = is1max; - // Ensure range: 1 <= is2 <= is2max - if ( divresult.rem > 1 ) - is2 = divresult.rem; - else - is2 = 1; - if ( is2max < is2 ) - is2 = is2max; - - int i = ((is1/177) % 177) + 2; - int j = (is1 % 177) + 2; - int k = ((is2/169) % 178) + 1; - int l = (is2 % 169); - - for (int ii = 0; ii < 97; ii++) { - double s = 0.0; - double t = 0.5; - for (int jj = 0; jj < 24; jj++) { - //m = mod(mod(i*j, 179)*k, 179) - int mtemp = (i*j) % 179; - int m = (mtemp*k) % 179; - i = j; - j = k; - k = m; - l = (((53*l) + 1) % 169); - //if (mod(l*m, 64) .ge. 32) s = s + t - if ( ((l*m) % 64) >= 32 ) s += t; - //t = 0.5d0 * t - t *= 0.5; - } - RN_generator.u[ii] = s; +/** \return Description of current default RNG. */ +const char* Random_Number::CurrentDefaultRngStr() { + const char* str = 0; + switch (defaultType_) { + case MARSAGLIA : str = "Marsaglia"; break; + case STDLIB : str = "C stdlib"; break; + case MERSENNE_TWISTER : str = "Mersenne Twister (mt19937)"; break; + case PCG32 : str = "Permuted Congruential Generator (32 bit)"; break; + case XOSHIRO128PP : str = "Xoshiro128++"; break; } + return str; +} - RN_generator.c = 362436.0 / 16777216.0; - RN_generator.cd = 7654321.0 / 16777216.0; - RN_generator.cm = 16777213.0 / 16777216.0; +/** Allocate RNG according to the current default type. */ +void Random_Number::allocateRng() { + if (rng_ != 0) delete rng_; + rng_ = 0; + switch (defaultType_) { + case MARSAGLIA : + mprintf("\tRNG: Marsaglia\n"); + rng_ = new Cpptraj::RNG_Marsaglia(); + break; + case STDLIB : + mprintf("\tRNG: C stdlib\n"); + rng_ = new Cpptraj::RNG_Stdlib(); + break; + case MERSENNE_TWISTER : + mprintf("\tRNG: Mersenne twister\n"); + rng_ = new Cpptraj::RNG_MersenneTwister(); + break; + case PCG32 : + mprintf("\tRNG: PCG 32\n"); + rng_ = new Cpptraj::RNG_PCG32(); + break; + case XOSHIRO128PP : + mprintf("\tRNG: Xoshiro128++\n"); + rng_ = new Cpptraj::RNG_Xoshiro128pp(); + break; + } +} - RN_generator.i97 = 96; - RN_generator.j97 = 32; +/** Allocate RNG, initialize with given seed. */ +int Random_Number::rn_set(int seedIn) { + allocateRng(); + return rng_->Set_Seed( seedIn ); } -// Random_Number::rn_set() -/** When called with no seed, use Amber default. */ -void Random_Number::rn_set( ) { - rn_set( 71277 ); +/** Initialize with default seed. */ +//void Random_Number::rn_set() { +// rng_->Set_Seed(); +//} + +/** Generate random number between 0 and 1. */ +double Random_Number::rn_gen() const { + return rng_->Generate(); } -// Random_Number::rn_gen() -/** \author Portable Random number generator by George Marsaglia - * \author Amber 3.0 Rev A implementation by George Seibel - * - * This random number generator originally appeared in 'Toward a Universal - * Random Number Generator' by George Marsaglia and Arif Zaman. Florida - * State University Report: FSU-SCRI-87-50 (1987) - * - * It was later modified by F. James and published in 'A Review of Pseudo- - * random Number Generators' - * - * This is claimed to be the best known random number generator available. - * It passes ALL of the tests for random number generators and has a - * period of 2^144, is completely portable (gives bit identical results on - * all machines with at least 24-bit mantissas in the floating point - * representation). - * - * The algorithm is a combination of a Fibonacci sequence (with lags of 97 - * and 33, and operation "subtraction plus one, modulo one") and an - * "arithmetic sequence" (using subtraction). - * - * \return A random number between 0.0 and 1.0 - * \return -1.0 if the random number generator is not initialized. +/** Generate a random integer. */ +unsigned int Random_Number::rn_num() const { + return rng_->Number(); +} + +/** Generate a random integer on interval imin to imax. + * NOTE: Unlike rn_num(), this returns an integer to allow intervals of + * negative numbers. */ -double Random_Number::rn_gen() { - if (!IsSet()) { - mprinterr("Error: random number generator not initialized."); - return -1.0; +int Random_Number::rn_num_interval_signed(int imin, int imax) const { + int iwidth = imax - imin; + if (iwidth < 1) { + mprinterr("Internal Error: Random_Number::rn_num_interval_signed(): Interval width <= 0.\n"); + return 0; } - - double uni = RN_generator.u[RN_generator.i97] - RN_generator.u[RN_generator.j97]; - if (uni < 0.0) uni += 1.0; - RN_generator.u[RN_generator.i97] = uni; - RN_generator.i97--; - if (RN_generator.i97 == -1) RN_generator.i97 = 96; - RN_generator.j97--; - if (RN_generator.j97 == -1) RN_generator.j97 = 96; - RN_generator.c -= RN_generator.cd; - if (RN_generator.c < 0.0) RN_generator.c += RN_generator.cm; - uni -= RN_generator.c; - if (uni < 0.0) uni += 1.0; - return uni; + unsigned int uwidth = (unsigned int)iwidth + 1; + unsigned int umod = (rn_num() % uwidth); + return (int)umod + imin; } -// Random_Number::rn_gauss() -/** This is a version of amrand() that adds the constraint of a gaussian - * distribution, with mean "am" and standard deviation "sd". This - * requires rn_set() to have been called first, and "uses up" the - * same sequence that rn_gen() does. +/** Generate a random integer on interval umin to umax. */ +unsigned int Random_Number::rn_num_interval(unsigned int umin, unsigned int umax) const { + if (umax > umin) { + unsigned int uwidth = (umax - umin) + 1; + unsigned int umod = (rn_num() % uwidth); + return umod + umin; + } else { + mprinterr("Internal Error: Random_Number::rn_num_interval(): Interval max <= interval min.\n"); + } + return 0; +} + +/** This uses Generate() to generate random numbers between 0 and 1 + * in a Gaussian distribution, with mean "am" and standard deviation "sd". */ -double Random_Number::rn_gauss(double am, double sd) { +double Random_Number::GenerateGauss(double am, double sd) const { if (!IsSet()) { mprinterr("Error: random number generator not initialized."); return -1.0; @@ -162,10 +146,10 @@ double Random_Number::rn_gauss(double am, double sd) { // Get two random numbers, even on (-1,1): bool generate = true; while (generate) { - double uni = rn_gen(); + double uni = rng_->Generate(); double zeta1 = uni + uni - 1.0; - uni = rn_gen(); + uni = rng_->Generate(); double zeta2 = uni + uni - 1.0; double tmp1 = zeta1 * zeta1 + zeta2 * zeta2; @@ -177,3 +161,31 @@ double Random_Number::rn_gauss(double am, double sd) { } return 0.0; } + +/** Generate pseudo-random Gaussian sequence. */ // TODO deprecate this version +double Random_Number::rn_gauss(double am, double sd) const { + return GenerateGauss(am, sd); +} + +/** Use modern version of the Fisher-Yates shuffle to randomly reorder the + * given points. + */ +void Random_Number::ShufflePoints( std::vector& PointIndices ) const { + for (unsigned int i = PointIndices.size() - 1; i != 1; i--) + { // 0 <= j <= i + unsigned int j = rn_num_interval(0, i); + int temp = PointIndices[j]; + PointIndices[j] = PointIndices[i]; + PointIndices[i] = temp; + } +} + +/** \return true if RNG has been set up. */ +bool Random_Number::IsSet() const { + return rng_->IsSet(); +} + +/** \return Value of RNG seed. */ +int Random_Number::Seed() const { + return rng_->Seed(); +} diff --git a/src/Random.h b/src/Random.h index 226e0241bb..396323653b 100644 --- a/src/Random.h +++ b/src/Random.h @@ -1,61 +1,52 @@ #ifndef INC_RANDOM_H #define INC_RANDOM_H -// Class: Random_Number -/// Marsaglias random number generator as implemented in Amber 3.0 Rev A -/** \author George Seibel - * - * Adapted from fortran code in $AMBERHOME/src/pmemd/src/random.F90 - * - * This random number generator originally appeared in "Toward a Universal - * Random Number Generator" by George Marsaglia and Arif Zaman. Florida - * State University Report: FSU-SCRI-87-50 (1987) - * - * It was later modified by F. James and published in "A Review of Pseudo- - * random Number Generators" - * - * This is claimed to be the best known random number generator available. - * It passes ALL of the tests for random number generators and has a - * period of 2^144, is completely portable (gives bit identical results on - * all machines with at least 24-bit mantissas in the floating point - * representation). - * - * The algorithm is a combination of a Fibonacci sequence (with lags of 97 - * and 33, and operation "subtraction plus one, modulo one") and an - * "arithmetic sequence" (using subtraction). - */ +#include +namespace Cpptraj { + class RNG; +} +/// Wrapper around RNG class Random_Number { public: Random_Number(); - /// Initialize the random number generator with the given seed - void rn_set(int); + ~Random_Number(); + /// Different possible RNGs + enum RngType { MARSAGLIA = 0, STDLIB, MERSENNE_TWISTER, PCG32, XOSHIRO128PP }; + + static void SetDefaultRng(RngType); + + //static void SetDefaultSeed(int); + + static const char* CurrentDefaultRngStr(); + + /// Allocate and initialize the random number generator with the given seed + int rn_set(int); /// Initialize RN generator with 71277 (Amber default) - void rn_set(); + //void rn_set(); /// Generate a random number between 0.0 and 1.0 - double rn_gen(); - /// Generate a pseudo-random Gaussian sequence. - double rn_gauss(double,double); + double rn_gen() const; + /// Generate a random integer + unsigned int rn_num() const; + /// Generate a random number on defined interval + int rn_num_interval_signed(int, int) const; + /// Generate a random number on defined interval + unsigned int rn_num_interval(unsigned int, unsigned int) const; + /// Generate random numbers between 0 and 1 in Gaussian distribution. + double GenerateGauss(double, double) const; + /// Generate a pseudo-random Gaussian sequence. // TODO deprecate this version + double rn_gauss(double,double) const; + /// Shuffle given points in array + void ShufflePoints(std::vector&) const; /// \return true if RN generator has been set up. - bool IsSet() const { return RN_generator.iseed != -1; } + bool IsSet() const; /// \return Value of RN generator seed - int Seed() const { return RN_generator.iseed; } + int Seed() const; private: - /// Variables necessary for Marsaglia random number stream. - /** This is placed in a struct in case the state of the random number - * generator ever needs to be stored. - */ - struct random_state { - // Real variables in Marsaglia algorithm - double u[97]; - double c; - double cd; - double cm; - // Pointers into u() in Marsaglia algorithm - int i97; - int j97; - // Random seed; if -1 random generator has not been set - int iseed; - }; - /// Hold the state of the random number generator - random_state RN_generator; + /// Allocate the RNG + void allocateRng(); + + Cpptraj::RNG* rng_; ///< Hold the random number generator implemenatation. + + static RngType defaultType_; ///< Default RNG type + //static int defaultSeed_; ///< Default RNG seed }; #endif diff --git a/src/TextFormat.cpp b/src/TextFormat.cpp index 2c97524e49..706c2d3e85 100644 --- a/src/TextFormat.cpp +++ b/src/TextFormat.cpp @@ -1,10 +1,10 @@ #include "TextFormat.h" #include "StringRoutines.h" -char TextFormat::TypeChar_[] = { 'f', 'E', 'g', 'i', 's' }; +char TextFormat::TypeChar_[] = { 'f', 'E', 'g', 'i', 'u', 's' }; const char* TextFormat::TypeDesc_[] = { - "DOUBLE", "SCIENTIFIC", "GENERAL", "INTEGER", "STRING" + "DOUBLE", "SCIENTIFIC", "GENERAL", "INTEGER", "UNSIGNED", "STRING" }; // TODO benchmark - will using a big buffer and C string routines be better? @@ -73,7 +73,7 @@ void TextFormat::SetCoordFormat(size_t maxFrames, double min, double step, // Default width for column is at least default_width. if (col_width < default_width) col_width = default_width; // Set column data format string, left-aligned (no leading space). - if (type_ == INTEGER || type_ == STRING) // sanity check + if (type_ == INTEGER || type_ == UNSIGNED || type_ == STRING) // sanity check type_ = DOUBLE; width_ = col_width; precision_ = col_precision; diff --git a/src/TextFormat.h b/src/TextFormat.h index 7998b3624a..dbf76f40b8 100644 --- a/src/TextFormat.h +++ b/src/TextFormat.h @@ -6,7 +6,7 @@ class TextFormat { public: /// Formats: f, E, g, i, s - enum FmtType { DOUBLE = 0, SCIENTIFIC, GDOUBLE, INTEGER, STRING }; + enum FmtType { DOUBLE = 0, SCIENTIFIC, GDOUBLE, INTEGER, UNSIGNED, STRING }; /// Alignment type: Right, Left, Right with leading space. enum AlignType { RIGHT = 0, LEFT, LEADING_SPACE }; /// CONSTRUCTOR - default 8.3f format diff --git a/src/Version.h b/src/Version.h index 1fe194b5c4..f1d309a9fe 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V5.0.5" +#define CPPTRAJ_INTERNAL_VERSION "V5.0.6" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 3bfd4c6b62..38e0c201a4 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -159,7 +159,7 @@ Cluster_ReadInfo.o : Cluster_ReadInfo.cpp ArgList.h ArrayIterator.h AssociatedDa Cmd.o : Cmd.cpp Cmd.h DispatchObject.h CmdInput.o : CmdInput.cpp CmdInput.h StringRoutines.h CmdList.o : CmdList.cpp Cmd.h CmdList.h DispatchObject.h -Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h ActionTopWriter.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_DihedralRMS.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_FixImagedBonds.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_HydrogenBond.h Action_Image.h Action_InfraredSpectrum.h Action_Jcoupling.h Action_LESsplit.h Action_LIE.h Action_LipidOrder.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Time.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Action_XtalSymm.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_Clustering.h Analysis_ConstantPHStats.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_EvalPlateau.h Analysis_FFT.h Analysis_HausdorffDistance.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Slope.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h AxisType.h BaseIOtype.h Box.h BoxArgs.h BufferedLine.h CharMask.h ClusterDist.h ClusterList.h ClusterMap.h ClusterNode.h ClusterSieve.h Cmd.h CmdInput.h CmdList.h Command.h CompactFrameArray.h ComplexArray.h Constants.h Constraints.h ControlBlock.h ControlBlock_For.h CoordinateInfo.h Corr.h Cph.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Deprecated.h DihedralSearch.h Dimension.h DispatchObject.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOutList.h ExclusionArray.h Exec.h Exec_AddMissingRes.h Exec_Analyze.h Exec_Calc.h Exec_CatCrd.h Exec_Change.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_CreateSet.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_Emin.h Exec_Flatten.h Exec_GenerateAmberRst.h Exec_Graft.h Exec_Help.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParallelAnalysis.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrepareForLeap.h Exec_PrintData.h Exec_ReadData.h Exec_ReadEnsembleData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_Set.h Exec_Show.h Exec_SortEnsembleData.h Exec_SplitCoords.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_UpdateParameters.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageOption.h ImageTypes.h InputTrajCommon.h MapAtom.h MaskArray.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h PubFFT.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Spline.h SplineFxnTable.h StructureCheck.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h cuda_kernels/GistCudaSetup.cuh molsurf.h +Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h ActionTopWriter.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_DihedralRMS.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_FixImagedBonds.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_HydrogenBond.h Action_Image.h Action_InfraredSpectrum.h Action_Jcoupling.h Action_LESsplit.h Action_LIE.h Action_LipidOrder.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Time.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Action_XtalSymm.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_Clustering.h Analysis_ConstantPHStats.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_EvalPlateau.h Analysis_FFT.h Analysis_HausdorffDistance.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Slope.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h AxisType.h BaseIOtype.h Box.h BoxArgs.h BufferedLine.h CharMask.h ClusterDist.h ClusterList.h ClusterMap.h ClusterNode.h ClusterSieve.h Cmd.h CmdInput.h CmdList.h Command.h CompactFrameArray.h ComplexArray.h Constants.h Constraints.h ControlBlock.h ControlBlock_For.h CoordinateInfo.h Corr.h Cph.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Deprecated.h DihedralSearch.h Dimension.h DispatchObject.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOutList.h ExclusionArray.h Exec.h Exec_AddMissingRes.h Exec_Analyze.h Exec_Calc.h Exec_CatCrd.h Exec_Change.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_CreateSet.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_Emin.h Exec_Flatten.h Exec_GenerateAmberRst.h Exec_Graft.h Exec_Help.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParallelAnalysis.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrepareForLeap.h Exec_PrintData.h Exec_Random.h Exec_ReadData.h Exec_ReadEnsembleData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_Set.h Exec_Show.h Exec_SortEnsembleData.h Exec_SplitCoords.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_UpdateParameters.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageOption.h ImageTypes.h InputTrajCommon.h MapAtom.h MaskArray.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h PubFFT.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Spline.h SplineFxnTable.h StructureCheck.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h cuda_kernels/GistCudaSetup.cuh molsurf.h CompactFrameArray.o : CompactFrameArray.cpp Box.h CompactFrameArray.h CoordinateInfo.h CpptrajStdio.h Matrix_3x3.h Parallel.h ReplicaDimArray.h Vec3.h ComplexArray.o : ComplexArray.cpp ArrayIterator.h ComplexArray.h Constraints.o : Constraints.cpp ArgList.h Atom.h AtomMask.h AtomType.h Box.h CharMask.h Constants.h Constraints.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -169,7 +169,7 @@ Corr.o : Corr.cpp ArrayIterator.h ComplexArray.h Corr.h PubFFT.h Cph.o : Cph.cpp Cph.h NameType.h Cpptraj.o : Cpptraj.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cmd.h CmdInput.h CmdList.h Command.h Constants.h CoordinateInfo.h Cpptraj.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ParmFile.h Range.h ReadLine.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h TopInfo.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Version.h CpptrajFile.o : CpptrajFile.cpp CpptrajFile.h CpptrajStdio.h FileIO.h FileIO_Bzip2.h FileIO_Gzip.h FileIO_Mpi.h FileIO_MpiShared.h FileIO_Std.h FileName.h Parallel.h StringRoutines.h -CpptrajState.o : CpptrajState.cpp Action.h ActionList.h ActionState.h Action_CreateCrd.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CompactFrameArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_Topology.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleNavigator.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h +CpptrajState.o : CpptrajState.cpp Action.h ActionList.h ActionState.h Action_CreateCrd.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CompactFrameArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_Topology.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleNavigator.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h CpptrajStdio.o : CpptrajStdio.cpp Parallel.h CurveFit.o : CurveFit.cpp CurveFit.h DataFile.o : DataFile.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h ClusterDist.h ClusterSieve.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataFile.h DataIO.h DataIO_CCP4.h DataIO_CharmmFastRep.h DataIO_CharmmOutput.h DataIO_CharmmRepLog.h DataIO_CharmmRtfPrm.h DataIO_Cmatrix.h DataIO_Cpout.h DataIO_Evecs.h DataIO_Gnuplot.h DataIO_Grace.h DataIO_Mdout.h DataIO_NC_Cmatrix.h DataIO_OpenDx.h DataIO_Peaks.h DataIO_RemLog.h DataIO_Std.h DataIO_VecTraj.h DataIO_XVG.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Cmatrix.h DataSet_Cmatrix_MEM.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_RemLog.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TypeNameHolder.h Unit.h Vec3.h @@ -195,7 +195,7 @@ DataIO_VecTraj.o : DataIO_VecTraj.cpp ActionFrameCounter.h ArgList.h ArrayIterat DataIO_XVG.o : DataIO_XVG.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_XVG.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataIO_Xplor.o : DataIO_Xplor.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataIO.h DataIO_Xplor.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h DataSet.o : DataSet.cpp AssociatedData.h CpptrajFile.h CpptrajStdio.h DataSet.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h -DataSetList.o : DataSetList.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h Box.h ClusterDist.h ClusterSieve.h CompactFrameArray.h ComplexArray.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_DISK.h DataSet_Cmatrix_MEM.h DataSet_Cmatrix_NOMEM.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_PHREMD.h DataSet_PHREMD_Explicit.h DataSet_PHREMD_Implicit.h DataSet_Parameters.h DataSet_RemLog.h DataSet_StringVar.h DataSet_Tensor.h DataSet_Topology.h DataSet_Vector.h DataSet_Vector_Scalar.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_disk.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Spline.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h SymmetricTensor.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Unit.h Vec3.h +DataSetList.o : DataSetList.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h Box.h ClusterDist.h ClusterSieve.h CompactFrameArray.h ComplexArray.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Cmatrix_DISK.h DataSet_Cmatrix_MEM.h DataSet_Cmatrix_NOMEM.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_PHREMD.h DataSet_PHREMD_Explicit.h DataSet_PHREMD_Implicit.h DataSet_Parameters.h DataSet_RemLog.h DataSet_StringVar.h DataSet_Tensor.h DataSet_Topology.h DataSet_Vector.h DataSet_Vector_Scalar.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_integer_disk.h DataSet_integer_mem.h DataSet_pH.h DataSet_string.h DataSet_unsignedInt.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h Hungarian.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NC_Cmatrix.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Spline.h StringRoutines.h SymbolExporting.h SymmetricRmsdCalc.h SymmetricTensor.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Unit.h Vec3.h DataSet_1D.o : DataSet_1D.cpp ArrayIterator.h AssociatedData.h ComplexArray.h Constants.h Corr.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h PubFFT.h Range.h TextFormat.h DataSet_3D.o : DataSet_3D.cpp AssociatedData.h Box.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_3D.h Dimension.h FileIO.h FileName.h GridBin.h Matrix_3x3.h MetaData.h Parallel.h Range.h TextFormat.h Vec3.h DataSet_Cmatrix.o : DataSet_Cmatrix.cpp ArrayIterator.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h Box.h ClusterDist.h ClusterSieve.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h DataSet_Cmatrix.h DataSet_Coords.h Dimension.h FileIO.h FileName.h Frame.h Hungarian.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h SymmetricRmsdCalc.h TextFormat.h Topology.h TypeNameHolder.h Unit.h Vec3.h @@ -228,6 +228,7 @@ DataSet_integer_disk.o : DataSet_integer_disk.cpp AssociatedData.h CpptrajFile.h DataSet_integer_mem.o : DataSet_integer_mem.cpp AssociatedData.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_integer.h DataSet_integer_mem.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h DataSet_pH.o : DataSet_pH.cpp AssociatedData.h Cph.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_1D.h DataSet_pH.h Dimension.h FileIO.h FileName.h MetaData.h NameType.h Parallel.h Range.h TextFormat.h DataSet_string.o : DataSet_string.cpp AssociatedData.h CpptrajFile.h DataSet.h DataSet_string.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h +DataSet_unsignedInt.o : DataSet_unsignedInt.cpp AssociatedData.h CpptrajFile.h DataSet.h DataSet_1D.h DataSet_unsignedInt.h Dimension.h FileIO.h FileName.h MetaData.h Parallel.h Range.h TextFormat.h Deprecated.o : Deprecated.cpp CpptrajStdio.h Deprecated.h DispatchObject.h DihedralSearch.o : DihedralSearch.cpp ArgList.h Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h DihedralSearch.h FileName.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h DistRoutines.o : DistRoutines.cpp Box.h DistRoutines.h ImageOption.h Matrix_3x3.h Parallel.h Vec3.h @@ -277,6 +278,7 @@ Exec_PermuteDihedrals.o : Exec_PermuteDihedrals.cpp Action.h ActionFrameCounter. Exec_Precision.o : Exec_Precision.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Precision.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_PrepareForLeap.o : Exec_PrepareForLeap.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_PrepareForLeap.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h ImageOption.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_PrintData.o : Exec_PrintData.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_PrintData.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h +Exec_Random.o : Exec_Random.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Random.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Random.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_ReadData.o : Exec_ReadData.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ReadData.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_ReadEnsembleData.o : Exec_ReadEnsembleData.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ReadEnsembleData.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_ReadInput.o : Exec_ReadInput.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cmd.h CmdList.h Command.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ReadInput.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h @@ -293,6 +295,7 @@ Exec_Top.o : Exec_Top.cpp Action.h ActionList.h ActionState.h Analysis.h Analysi Exec_Traj.o : Exec_Traj.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Traj.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_UpdateParameters.o : Exec_UpdateParameters.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Parameters.h DataSet_Topology.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_UpdateParameters.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h Exec_ViewRst.o : Exec_ViewRst.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h ViewRst.h +ExternalFxn.o : ExternalFxn.cpp Random.h FileIO_Bzip2.o : FileIO_Bzip2.cpp CpptrajStdio.h FileIO.h FileIO_Bzip2.h FileIO_Gzip.o : FileIO_Gzip.cpp CpptrajStdio.h FileIO.h FileIO_Gzip.h FileIO_Mpi.o : FileIO_Mpi.cpp FileIO.h FileIO_Mpi.h Parallel.h @@ -357,8 +360,14 @@ PotentialTerm_OpenMM.o : PotentialTerm_OpenMM.cpp Atom.h AtomMask.h AtomType.h B ProgressBar.o : ProgressBar.cpp CpptrajStdio.h ProgressBar.h ProgressTimer.o : ProgressTimer.cpp CpptrajStdio.h ProgressTimer.h Timer.h PubFFT.o : PubFFT.cpp ArrayIterator.h ComplexArray.h CpptrajStdio.h PubFFT.h +RNG.o : RNG.cpp CpptrajStdio.h RNG.h +RNG_Marsaglia.o : RNG_Marsaglia.cpp CpptrajStdio.h RNG.h RNG_Marsaglia.h +RNG_MersenneTwister.o : RNG_MersenneTwister.cpp CpptrajStdio.h RNG.h RNG_MersenneTwister.h +RNG_PCG32.o : RNG_PCG32.cpp CpptrajStdio.h RNG.h RNG_PCG32.h pcg_random.hpp +RNG_Stdlib.o : RNG_Stdlib.cpp RNG.h RNG_Stdlib.h +RNG_Xoshiro128pp.o : RNG_Xoshiro128pp.cpp RNG.h RNG_Xoshiro128pp.h xoshiro128plusplus.h RPNcalc.o : RPNcalc.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_MatrixDbl.h DataSet_Mesh.h DataSet_Vector.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h -Random.o : Random.cpp CpptrajStdio.h Random.h +Random.o : Random.cpp CpptrajStdio.h RNG.h RNG_Marsaglia.h RNG_MersenneTwister.h RNG_PCG32.h RNG_Stdlib.h RNG_Xoshiro128pp.h Random.h pcg_random.hpp Range.o : Range.cpp ArgList.h CpptrajStdio.h Range.h ReadLine.o : ReadLine.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cmd.h CmdInput.h CmdList.h Command.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReadLine.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h ReferenceAction.o : ReferenceAction.cpp ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Coords_TRJ.h Dimension.h FileIO.h FileName.h Frame.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajFrameIndex.h Trajin.h TypeNameHolder.h Unit.h Vec3.h @@ -416,3 +425,4 @@ ViewRst.o : ViewRst.cpp ActionFrameCounter.h ArgList.h AssociatedData.h Atom.h A main.o : main.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h Cpptraj.h CpptrajFile.h CpptrajState.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h molsurf.o : molsurf.c molsurf.h vmdplugin/dtrplugin.o : vmdplugin/dtrplugin.cpp vmdplugin/../ByteRoutines.h vmdplugin/dtrplugin.hxx vmdplugin/vmddir.h +xoshiro128plusplus.o : xoshiro128plusplus.cpp diff --git a/src/cpptrajfiles b/src/cpptrajfiles index 2047c673e3..f4adab6365 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -229,6 +229,7 @@ COMMON_SOURCES= \ DataSet_string.cpp \ DataSet_StringVar.cpp \ DataSet_Tensor.cpp \ + DataSet_unsignedInt.cpp \ Deprecated.cpp \ DihedralSearch.cpp \ DistRoutines.cpp \ @@ -277,6 +278,7 @@ COMMON_SOURCES= \ Exec_Precision.cpp \ Exec_PrepareForLeap.cpp \ Exec_PrintData.cpp \ + Exec_Random.cpp \ Exec_ReadData.cpp \ Exec_ReadEnsembleData.cpp \ Exec_ReadInput.cpp \ @@ -357,8 +359,15 @@ COMMON_SOURCES= \ ProgressBar.cpp \ ProgressTimer.cpp \ PubFFT.cpp \ + ExternalFxn.cpp \ Random.cpp \ Range.cpp \ + RNG.cpp \ + RNG_Marsaglia.cpp \ + RNG_MersenneTwister.cpp \ + RNG_PCG32.cpp \ + RNG_Stdlib.cpp \ + RNG_Xoshiro128pp.cpp \ RPNcalc.cpp \ ReferenceAction.cpp \ RemdReservoirNC.cpp \ @@ -411,9 +420,10 @@ COMMON_SOURCES= \ TrajoutList.cpp \ Vec3.cpp \ ViewRst.cpp \ - vmdplugin/dtrplugin.cpp + vmdplugin/dtrplugin.cpp \ + xoshiro128plusplus.cpp -CSOURCES= molsurf.c +CSOURCES= molsurf.c # The files below are used by cpptraj and libcpptraj but must be compiled # differently for libcpptraj. diff --git a/src/pcg_extras.hpp b/src/pcg_extras.hpp new file mode 100644 index 0000000000..ec3e5694e2 --- /dev/null +++ b/src/pcg_extras.hpp @@ -0,0 +1,637 @@ +/* + * PCG Random Number Generation for C++ + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This file provides support code that is useful for random-number generation + * but not specific to the PCG generation scheme, including: + * - 128-bit int support for platforms where it isn't available natively + * - bit twiddling operations + * - I/O of 128-bit and 8-bit integers + * - Handling the evilness of SeedSeq + * - Support for efficiently producing random numbers less than a given + * bound + */ + +#ifndef PCG_EXTRAS_HPP_INCLUDED +#define PCG_EXTRAS_HPP_INCLUDED 1 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef __GNUC__ + #include +#endif + +/* + * Abstractions for compiler-specific directives + */ + +#ifdef __GNUC__ + #define PCG_NOINLINE __attribute__((noinline)) +#else + #define PCG_NOINLINE +#endif + +/* + * Some members of the PCG library use 128-bit math. When compiling on 64-bit + * platforms, both GCC and Clang provide 128-bit integer types that are ideal + * for the job. + * + * On 32-bit platforms (or with other compilers), we fall back to a C++ + * class that provides 128-bit unsigned integers instead. It may seem + * like we're reinventing the wheel here, because libraries already exist + * that support large integers, but most existing libraries provide a very + * generic multiprecision code, but here we're operating at a fixed size. + * Also, most other libraries are fairly heavyweight. So we use a direct + * implementation. Sadly, it's much slower than hand-coded assembly or + * direct CPU support. + * + */ +#if __SIZEOF_INT128__ + namespace pcg_extras { + typedef __uint128_t pcg128_t; + } + #define PCG_128BIT_CONSTANT(high,low) \ + ((pcg128_t(high) << 64) + low) +#else + #include "pcg_uint128.hpp" + namespace pcg_extras { + typedef pcg_extras::uint_x4 pcg128_t; + } + #define PCG_128BIT_CONSTANT(high,low) \ + pcg128_t(high,low) + #define PCG_EMULATED_128BIT_MATH 1 +#endif + + +namespace pcg_extras { + +/* + * We often need to represent a "number of bits". When used normally, these + * numbers are never greater than 128, so an unsigned char is plenty. + * If you're using a nonstandard generator of a larger size, you can set + * PCG_BITCOUNT_T to have it define it as a larger size. (Some compilers + * might produce faster code if you set it to an unsigned int.) + */ + +#ifndef PCG_BITCOUNT_T + typedef uint8_t bitcount_t; +#else + typedef PCG_BITCOUNT_T bitcount_t; +#endif + +/* + * C++ requires us to be able to serialize RNG state by printing or reading + * it from a stream. Because we use 128-bit ints, we also need to be able + * ot print them, so here is code to do so. + * + * This code provides enough functionality to print 128-bit ints in decimal + * and zero-padded in hex. It's not a full-featured implementation. + */ + +template +std::basic_ostream& +operator<<(std::basic_ostream& out, pcg128_t value) +{ + auto desired_base = out.flags() & out.basefield; + bool want_hex = desired_base == out.hex; + + if (want_hex) { + uint64_t highpart = uint64_t(value >> 64); + uint64_t lowpart = uint64_t(value); + auto desired_width = out.width(); + if (desired_width > 16) { + out.width(desired_width - 16); + } + if (highpart != 0 || desired_width > 16) + out << highpart; + CharT oldfill; + if (highpart != 0) { + out.width(16); + oldfill = out.fill('0'); + } + auto oldflags = out.setf(decltype(desired_base){}, out.showbase); + out << lowpart; + out.setf(oldflags); + if (highpart != 0) { + out.fill(oldfill); + } + return out; + } + constexpr size_t MAX_CHARS_128BIT = 40; + + char buffer[MAX_CHARS_128BIT]; + char* pos = buffer+sizeof(buffer); + *(--pos) = '\0'; + constexpr auto BASE = pcg128_t(10ULL); + do { + auto div = value / BASE; + auto mod = uint32_t(value - (div * BASE)); + *(--pos) = '0' + mod; + value = div; + } while(value != pcg128_t(0ULL)); + return out << pos; +} + +template +std::basic_istream& +operator>>(std::basic_istream& in, pcg128_t& value) +{ + typename std::basic_istream::sentry s(in); + + if (!s) + return in; + + constexpr auto BASE = pcg128_t(10ULL); + pcg128_t current(0ULL); + bool did_nothing = true; + bool overflow = false; + for(;;) { + CharT wide_ch = in.get(); + if (!in.good()) + break; + auto ch = in.narrow(wide_ch, '\0'); + if (ch < '0' || ch > '9') { + in.unget(); + break; + } + did_nothing = false; + pcg128_t digit(uint32_t(ch - '0')); + pcg128_t timesbase = current*BASE; + overflow = overflow || timesbase < current; + current = timesbase + digit; + overflow = overflow || current < digit; + } + + if (did_nothing || overflow) { + in.setstate(std::ios::failbit); + if (overflow) + current = ~pcg128_t(0ULL); + } + + value = current; + + return in; +} + +/* + * Likewise, if people use tiny rngs, we'll be serializing uint8_t. + * If we just used the provided IO operators, they'd read/write chars, + * not ints, so we need to define our own. We *can* redefine this operator + * here because we're in our own namespace. + */ + +template +std::basic_ostream& +operator<<(std::basic_ostream&out, uint8_t value) +{ + return out << uint32_t(value); +} + +template +std::basic_istream& +operator>>(std::basic_istream& in, uint8_t target) +{ + uint32_t value = 0xdecea5edU; + in >> value; + if (!in && value == 0xdecea5edU) + return in; + if (value > uint8_t(~0)) { + in.setstate(std::ios::failbit); + value = ~0U; + } + target = uint8_t(value); + return in; +} + +/* Unfortunately, the above functions don't get found in preference to the + * built in ones, so we create some more specific overloads that will. + * Ugh. + */ + +inline std::ostream& operator<<(std::ostream& out, uint8_t value) +{ + return pcg_extras::operator<< (out, value); +} + +inline std::istream& operator>>(std::istream& in, uint8_t& value) +{ + return pcg_extras::operator>> (in, value); +} + + + +/* + * Useful bitwise operations. + */ + +/* + * XorShifts are invertable, but they are someting of a pain to invert. + * This function backs them out. It's used by the whacky "inside out" + * generator defined later. + */ + +template +inline itype unxorshift(itype x, bitcount_t bits, bitcount_t shift) +{ + if (2*shift >= bits) { + return x ^ (x >> shift); + } + itype lowmask1 = (itype(1U) << (bits - shift*2)) - 1; + itype highmask1 = ~lowmask1; + itype top1 = x; + itype bottom1 = x & lowmask1; + top1 ^= top1 >> shift; + top1 &= highmask1; + x = top1 | bottom1; + itype lowmask2 = (itype(1U) << (bits - shift)) - 1; + itype bottom2 = x & lowmask2; + bottom2 = unxorshift(bottom2, bits - shift, shift); + bottom2 &= lowmask1; + return top1 | bottom2; +} + +/* + * Rotate left and right. + * + * In ideal world, compilers would spot idiomatic rotate code and convert it + * to a rotate instruction. Of course, opinions vary on what the correct + * idiom is and how to spot it. For clang, sometimes it generates better + * (but still crappy) code if you define PCG_USE_ZEROCHECK_ROTATE_IDIOM. + */ + +template +inline itype rotl(itype value, bitcount_t rot) +{ + constexpr bitcount_t bits = sizeof(itype) * 8; + constexpr bitcount_t mask = bits - 1; +#if PCG_USE_ZEROCHECK_ROTATE_IDIOM + return rot ? (value << rot) | (value >> (bits - rot)) : value; +#else + return (value << rot) | (value >> ((- rot) & mask)); +#endif +} + +template +inline itype rotr(itype value, bitcount_t rot) +{ + constexpr bitcount_t bits = sizeof(itype) * 8; + constexpr bitcount_t mask = bits - 1; +#if PCG_USE_ZEROCHECK_ROTATE_IDIOM + return rot ? (value >> rot) | (value << (bits - rot)) : value; +#else + return (value >> rot) | (value << ((- rot) & mask)); +#endif +} + +/* Unfortunately, both Clang and GCC sometimes perform poorly when it comes + * to properly recognizing idiomatic rotate code, so for we also provide + * assembler directives (enabled with PCG_USE_INLINE_ASM). Boo, hiss. + * (I hope that these compilers get better so that this code can die.) + * + * These overloads will be preferred over the general template code above. + */ +#if PCG_USE_INLINE_ASM && __GNUC__ && (__x86_64__ || __i386__) + +inline uint8_t rotr(uint8_t value, bitcount_t rot) +{ + asm ("rorb %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +} + +inline uint16_t rotr(uint16_t value, bitcount_t rot) +{ + asm ("rorw %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +} + +inline uint32_t rotr(uint32_t value, bitcount_t rot) +{ + asm ("rorl %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +} + +#if __x86_64__ +inline uint64_t rotr(uint64_t value, bitcount_t rot) +{ + asm ("rorq %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +} +#endif // __x86_64__ + +#endif // PCG_USE_INLINE_ASM + + +/* + * The C++ SeedSeq concept (modelled by seed_seq) can fill an array of + * 32-bit integers with seed data, but sometimes we want to produce + * larger or smaller integers. + * + * The following code handles this annoyance. + * + * uneven_copy will copy an array of 32-bit ints to an array of larger or + * smaller ints (actually, the code is general it only needing forward + * iterators). The copy is identical to the one that would be performed if + * we just did memcpy on a standard little-endian machine, but works + * regardless of the endian of the machine (or the weirdness of the ints + * involved). + * + * generate_to initializes an array of integers using a SeedSeq + * object. It is given the size as a static constant at compile time and + * tries to avoid memory allocation. If we're filling in 32-bit constants + * we just do it directly. If we need a separate buffer and it's small, + * we allocate it on the stack. Otherwise, we fall back to heap allocation. + * Ugh. + * + * generate_one produces a single value of some integral type using a + * SeedSeq object. + */ + + /* uneven_copy helper, case where destination ints are less than 32 bit. */ + +template +SrcIter uneven_copy_impl( + SrcIter src_first, DestIter dest_first, DestIter dest_last, + std::true_type) +{ + typedef typename std::iterator_traits::value_type src_t; + typedef typename std::iterator_traits::value_type dest_t; + + constexpr bitcount_t SRC_SIZE = sizeof(src_t); + constexpr bitcount_t DEST_SIZE = sizeof(dest_t); + constexpr bitcount_t DEST_BITS = DEST_SIZE * 8; + constexpr bitcount_t SCALE = SRC_SIZE / DEST_SIZE; + + size_t count = 0; + src_t value; + + while (dest_first != dest_last) { + if ((count++ % SCALE) == 0) + value = *src_first++; // Get more bits + else + value >>= DEST_BITS; // Move down bits + + *dest_first++ = dest_t(value); // Truncates, ignores high bits. + } + return src_first; +} + + /* uneven_copy helper, case where destination ints are more than 32 bit. */ + +template +SrcIter uneven_copy_impl( + SrcIter src_first, DestIter dest_first, DestIter dest_last, + std::false_type) +{ + typedef typename std::iterator_traits::value_type src_t; + typedef typename std::iterator_traits::value_type dest_t; + + constexpr auto SRC_SIZE = sizeof(src_t); + constexpr auto SRC_BITS = SRC_SIZE * 8; + constexpr auto DEST_SIZE = sizeof(dest_t); + constexpr auto SCALE = (DEST_SIZE+SRC_SIZE-1) / SRC_SIZE; + + while (dest_first != dest_last) { + dest_t value(0UL); + unsigned int shift = 0; + + for (size_t i = 0; i < SCALE; ++i) { + value |= dest_t(*src_first++) << shift; + shift += SRC_BITS; + } + + *dest_first++ = value; + } + return src_first; +} + +/* uneven_copy, call the right code for larger vs. smaller */ + +template +inline SrcIter uneven_copy(SrcIter src_first, + DestIter dest_first, DestIter dest_last) +{ + typedef typename std::iterator_traits::value_type src_t; + typedef typename std::iterator_traits::value_type dest_t; + + constexpr bool DEST_IS_SMALLER = sizeof(dest_t) < sizeof(src_t); + + return uneven_copy_impl(src_first, dest_first, dest_last, + std::integral_constant{}); +} + +/* generate_to, fill in a fixed-size array of integral type using a SeedSeq + * (actually works for any random-access iterator) + */ + +template +inline void generate_to_impl(SeedSeq&& generator, DestIter dest, + std::true_type) +{ + generator.generate(dest, dest+size); +} + +template +void generate_to_impl(SeedSeq&& generator, DestIter dest, + std::false_type) +{ + typedef typename std::iterator_traits::value_type dest_t; + constexpr auto DEST_SIZE = sizeof(dest_t); + constexpr auto GEN_SIZE = sizeof(uint32_t); + + constexpr bool GEN_IS_SMALLER = GEN_SIZE < DEST_SIZE; + constexpr size_t FROM_ELEMS = + GEN_IS_SMALLER + ? size * ((DEST_SIZE+GEN_SIZE-1) / GEN_SIZE) + : (size + (GEN_SIZE / DEST_SIZE) - 1) + / ((GEN_SIZE / DEST_SIZE) + GEN_IS_SMALLER); + // this odd code ^^^^^^^^^^^^^^^^^ is work-around for + // a bug: http://llvm.org/bugs/show_bug.cgi?id=21287 + + if (FROM_ELEMS <= 1024) { + uint32_t buffer[FROM_ELEMS]; + generator.generate(buffer, buffer+FROM_ELEMS); + uneven_copy(buffer, dest, dest+size); + } else { + uint32_t* buffer = (uint32_t*) malloc(GEN_SIZE * FROM_ELEMS); + generator.generate(buffer, buffer+FROM_ELEMS); + uneven_copy(buffer, dest, dest+size); + free(buffer); + } +} + +template +inline void generate_to(SeedSeq&& generator, DestIter dest) +{ + typedef typename std::iterator_traits::value_type dest_t; + constexpr bool IS_32BIT = sizeof(dest_t) == sizeof(uint32_t); + + generate_to_impl(std::forward(generator), dest, + std::integral_constant{}); +} + +/* generate_one, produce a value of integral type using a SeedSeq + * (optionally, we can have it produce more than one and pick which one + * we want) + */ + +template +inline UInt generate_one(SeedSeq&& generator) +{ + UInt result[N]; + generate_to(std::forward(generator), result); + return result[i]; +} + +template +auto bounded_rand(RngType& rng, typename RngType::result_type upper_bound) + -> typename RngType::result_type +{ + typedef typename RngType::result_type rtype; + rtype threshold = (RngType::max() - RngType::min() + rtype(1) - upper_bound) + % upper_bound; + for (;;) { + rtype r = rng() - RngType::min(); + if (r >= threshold) + return r % upper_bound; + } +} + +template +void shuffle(Iter from, Iter to, RandType&& rng) +{ + typedef typename std::iterator_traits::difference_type delta_t; + auto count = to - from; + while (count > 1) { + delta_t chosen(bounded_rand(rng, count)); + --count; + --to; + using std::swap; + swap(*(from+chosen), *to); + } +} + +/* + * Although std::seed_seq is useful, it isn't everything. Often we want to + * initialize a random-number generator some other way, such as from a random + * device. + * + * Technically, it does not meet the requirements of a SeedSequence because + * it lacks some of the rarely-used member functions (some of which would + * be impossible to provide). However the C++ standard is quite specific + * that actual engines only called the generate method, so it ought not to be + * a problem in practice. + */ + +template +class seed_seq_from { +private: + RngType rng_; + + typedef uint_least32_t result_type; + +public: + template + seed_seq_from(Args&&... args) : + rng_(std::forward(args)...) + { + // Nothing (else) to do... + } + + template + void generate(Iter start, Iter finish) + { + for (auto i = start; i != finish; ++i) + *i = result_type(rng_()); + } + + constexpr size_t size() const + { + return (sizeof(typename RngType::result_type) > sizeof(result_type) + && RngType::max() > ~size_t(0UL)) + ? ~size_t(0UL) + : size_t(RngType::max()); + } +}; + +/* + * Sometimes you might want a distinct seed based on when the program + * was compiled. That way, a particular instance of the program will + * behave the same way, but when recompiled it'll produce a different + * value. + */ + +template +struct static_arbitrary_seed { +private: + static constexpr IntType fnv(IntType hash, const char* pos) { + return *pos == '\0' + ? hash + : fnv((hash * IntType(16777619U)) ^ *pos, (pos+1)); + } + +public: + static constexpr IntType value = fnv(IntType(2166136261U ^ sizeof(IntType)), + __DATE__ __TIME__ __FILE__); +}; + +// Sometimes, when debugging or testing, it's handy to be able print the name +// of a (in human-readable form). This code allows the idiom: +// +// cout << printable_typename() +// +// to print out my_foo_type_t (or its concrete type if it is a synonym) + +template +struct printable_typename {}; + +template +std::ostream& operator<<(std::ostream& out, printable_typename) { + const char *implementation_typename = typeid(T).name(); +#ifdef __GNUC__ + int status; + const char* pretty_name = + abi::__cxa_demangle(implementation_typename, NULL, NULL, &status); + if (status == 0) + out << pretty_name; + free((void*) pretty_name); + if (status == 0) + return out; +#endif + out << implementation_typename; + return out; +} + +} // namespace pcg_extras + +#endif // PCG_EXTRAS_HPP_INCLUDED diff --git a/src/pcg_random.hpp b/src/pcg_random.hpp new file mode 100644 index 0000000000..3f04d854e6 --- /dev/null +++ b/src/pcg_random.hpp @@ -0,0 +1,1751 @@ +/* + * PCG Random Number Generation for C++ + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code provides the reference implementation of the PCG family of + * random number generators. The code is complex because it implements + * + * - several members of the PCG family, specifically members corresponding + * to the output functions: + * - XSH RR (good for 64-bit state, 32-bit output) + * - XSH RS (good for 64-bit state, 32-bit output) + * - XSL RR (good for 128-bit state, 64-bit output) + * - RXS M XS (statistically most powerful generator) + * - XSL RR RR (good for 128-bit state, 128-bit output) + * - and RXS, RXS M, XSH, XSL (mostly for testing) + * - at potentially *arbitrary* bit sizes + * - with four different techniques for random streams (MCG, one-stream + * LCG, settable-stream LCG, unique-stream LCG) + * - and the extended generation schemes allowing arbitrary periods + * - with all features of C++11 random number generation (and more), + * some of which are somewhat painful, including + * - initializing with a SeedSequence which writes 32-bit values + * to memory, even though the state of the generator may not + * use 32-bit values (it might use smaller or larger integers) + * - I/O for RNGs and a prescribed format, which needs to handle + * the issue that 8-bit and 128-bit integers don't have working + * I/O routines (e.g., normally 8-bit = char, not integer) + * - equality and inequality for RNGs + * - and a number of convenience typedefs to mask all the complexity + * + * The code employes a fairly heavy level of abstraction, and has to deal + * with various C++ minutia. If you're looking to learn about how the PCG + * scheme works, you're probably best of starting with one of the other + * codebases (see www.pcg-random.org). But if you're curious about the + * constants for the various output functions used in those other, simpler, + * codebases, this code shows how they are calculated. + * + * On the positive side, at least there are convenience typedefs so that you + * can say + * + * pcg32 myRNG; + * + * rather than: + * + * pcg_detail::engine< + * uint32_t, // Output Type + * uint64_t, // State Type + * pcg_detail::xsh_rr_mixin, true, // Output Func + * pcg_detail::specific_stream, // Stream Kind + * pcg_detail::default_multiplier // LCG Mult + * > myRNG; + * + */ + +#ifndef PCG_RAND_HPP_INCLUDED +#define PCG_RAND_HPP_INCLUDED 1 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* + * The pcg_extras namespace contains some support code that is likley to + * be useful for a variety of RNGs, including: + * - 128-bit int support for platforms where it isn't available natively + * - bit twiddling operations + * - I/O of 128-bit and 8-bit integers + * - Handling the evilness of SeedSeq + * - Support for efficiently producing random numbers less than a given + * bound + */ + +#include "pcg_extras.hpp" + +namespace pcg_detail { + +using namespace pcg_extras; + +/* + * The LCG generators need some constants to function. This code lets you + * look up the constant by *type*. For example + * + * default_multiplier::multiplier() + * + * gives you the default multipler for 32-bit integers. We use the name + * of the constant and not a generic word like value to allow these classes + * to be used as mixins. + */ + +template +struct default_multiplier { + // Not defined for an arbitrary type +}; + +template +struct default_increment { + // Not defined for an arbitrary type +}; + +#define PCG_DEFINE_CONSTANT(type, what, kind, constant) \ + template <> \ + struct what ## _ ## kind { \ + static constexpr type kind() { \ + return constant; \ + } \ + }; + +PCG_DEFINE_CONSTANT(uint8_t, default, multiplier, 141U) +PCG_DEFINE_CONSTANT(uint8_t, default, increment, 77U) + +PCG_DEFINE_CONSTANT(uint16_t, default, multiplier, 12829U) +PCG_DEFINE_CONSTANT(uint16_t, default, increment, 47989U) + +PCG_DEFINE_CONSTANT(uint32_t, default, multiplier, 747796405U) +PCG_DEFINE_CONSTANT(uint32_t, default, increment, 2891336453U) + +PCG_DEFINE_CONSTANT(uint64_t, default, multiplier, 6364136223846793005ULL) +PCG_DEFINE_CONSTANT(uint64_t, default, increment, 1442695040888963407ULL) + +PCG_DEFINE_CONSTANT(pcg128_t, default, multiplier, + PCG_128BIT_CONSTANT(2549297995355413924ULL,4865540595714422341ULL)) +PCG_DEFINE_CONSTANT(pcg128_t, default, increment, + PCG_128BIT_CONSTANT(6364136223846793005ULL,1442695040888963407ULL)) + + +/* + * Each PCG generator is available in four variants, based on how it applies + * the additive constant for its underlying LCG; the variations are: + * + * single stream - all instances use the same fixed constant, thus + * the RNG always somewhere in same sequence + * mcg - adds zero, resulting in a single stream and reduced + * period + * specific stream - the constant can be changed at any time, selecting + * a different random sequence + * unique stream - the constant is based on the memory addresss of the + * object, thus every RNG has its own unique sequence + * + * This variation is provided though mixin classes which define a function + * value called increment() that returns the nesessary additive constant. + */ + + + +/* + * unique stream + */ + + +template +class unique_stream { +protected: + static constexpr bool is_mcg = false; + + // Is never called, but is provided for symmetry with specific_stream + void set_stream(...) + { + abort(); + } + +public: + typedef itype state_type; + + constexpr itype increment() const { + return itype(reinterpret_cast(this) | 1); + } + + constexpr itype stream() const + { + return increment() >> 1; + } + + static constexpr bool can_specify_stream = false; + + static constexpr size_t streams_pow2() + { + return (sizeof(itype) < sizeof(size_t) ? sizeof(itype) + : sizeof(size_t))*8 - 1u; + } + +protected: + constexpr unique_stream() = default; +}; + + +/* + * no stream (mcg) + */ + +template +class no_stream { +protected: + static constexpr bool is_mcg = true; + + // Is never called, but is provided for symmetry with specific_stream + void set_stream(...) + { + abort(); + } + +public: + typedef itype state_type; + + static constexpr itype increment() { + return 0; + } + + static constexpr bool can_specify_stream = false; + + static constexpr size_t streams_pow2() + { + return 0u; + } + +protected: + constexpr no_stream() = default; +}; + + +/* + * single stream/sequence (oneseq) + */ + +template +class oneseq_stream : public default_increment { +protected: + static constexpr bool is_mcg = false; + + // Is never called, but is provided for symmetry with specific_stream + void set_stream(...) + { + abort(); + } + +public: + typedef itype state_type; + + static constexpr itype stream() + { + return default_increment::increment() >> 1; + } + + static constexpr bool can_specify_stream = false; + + static constexpr size_t streams_pow2() + { + return 0u; + } + +protected: + constexpr oneseq_stream() = default; +}; + + +/* + * specific stream + */ + +template +class specific_stream { +protected: + static constexpr bool is_mcg = false; + + itype inc_ = default_increment::increment(); + +public: + typedef itype state_type; + typedef itype stream_state; + + constexpr itype increment() const { + return inc_; + } + + itype stream() + { + return inc_ >> 1; + } + + void set_stream(itype specific_seq) + { + inc_ = (specific_seq << 1) | 1; + } + + static constexpr bool can_specify_stream = true; + + static constexpr size_t streams_pow2() + { + return (sizeof(itype)*8) - 1u; + } + +protected: + specific_stream() = default; + + specific_stream(itype specific_seq) + : inc_((specific_seq << 1) | itype(1U)) + { + // Nothing (else) to do. + } +}; + + +/* + * This is where it all comes together. This function joins together three + * mixin classes which define + * - the LCG additive constant (the stream) + * - the LCG multiplier + * - the output function + * in addition, we specify the type of the LCG state, and the result type, + * and whether to use the pre-advance version of the state for the output + * (increasing instruction-level parallelism) or the post-advance version + * (reducing register pressure). + * + * Given the high level of parameterization, the code has to use some + * template-metaprogramming tricks to handle some of the suble variations + * involved. + */ + +template , + typename multiplier_mixin = default_multiplier > +class engine : protected output_mixin, + public stream_mixin, + protected multiplier_mixin { +protected: + itype state_; + + struct can_specify_stream_tag {}; + struct no_specifiable_stream_tag {}; + + using stream_mixin::increment; + using multiplier_mixin::multiplier; + +public: + typedef xtype result_type; + typedef itype state_type; + + static constexpr size_t period_pow2() + { + return sizeof(state_type)*8 - 2*stream_mixin::is_mcg; + } + + // It would be nice to use std::numeric_limits for these, but + // we can't be sure that it'd be defined for the 128-bit types. + + static constexpr result_type min() + { + return result_type(0UL); + } + + static constexpr result_type max() + { + return ~result_type(0UL); + } + +protected: + itype bump(itype state) + { + return state * multiplier() + increment(); + } + + itype base_generate() + { + return state_ = bump(state_); + } + + itype base_generate0() + { + itype old_state = state_; + state_ = bump(state_); + return old_state; + } + +public: + result_type operator()() + { + if (output_previous) + return this->output(base_generate0()); + else + return this->output(base_generate()); + } + + result_type operator()(result_type upper_bound) + { + return bounded_rand(*this, upper_bound); + } + +protected: + static itype advance(itype state, itype delta, + itype cur_mult, itype cur_plus); + + static itype distance(itype cur_state, itype newstate, itype cur_mult, + itype cur_plus, itype mask = ~itype(0U)); + + itype distance(itype newstate, itype mask = ~itype(0U)) const + { + return distance(state_, newstate, multiplier(), increment(), mask); + } + +public: + void advance(itype delta) + { + state_ = advance(state_, delta, this->multiplier(), this->increment()); + } + + void backstep(itype delta) + { + advance(-delta); + } + + void discard(itype delta) + { + advance(delta); + } + + bool wrapped() + { + if (stream_mixin::is_mcg) { + // For MCGs, the low order two bits never change. In this + // implementation, we keep them fixed at 3 to make this test + // easier. + return state_ == 3; + } else { + return state_ == 0; + } + } + + engine(itype state = itype(0xcafef00dd15ea5e5ULL)) + : state_(this->is_mcg ? state|state_type(3U) + : bump(state + this->increment())) + { + // Nothing else to do. + } + + // This function may or may not exist. It thus has to be a template + // to use SFINAE; users don't have to worry about its template-ness. + + template + engine(itype state, typename sm::stream_state stream_seed) + : stream_mixin(stream_seed), + state_(this->is_mcg ? state|state_type(3U) + : bump(state + this->increment())) + { + // Nothing else to do. + } + + template + engine(SeedSeq&& seedSeq, typename std::enable_if< + !stream_mixin::can_specify_stream + && !std::is_convertible::value + && !std::is_convertible::value, + no_specifiable_stream_tag>::type = {}) + : engine(generate_one(std::forward(seedSeq))) + { + // Nothing else to do. + } + + template + engine(SeedSeq&& seedSeq, typename std::enable_if< + stream_mixin::can_specify_stream + && !std::is_convertible::value + && !std::is_convertible::value, + can_specify_stream_tag>::type = {}) + : engine(generate_one(seedSeq), + generate_one(seedSeq)) + { + // Nothing else to do. + } + + + template + void seed(Args&&... args) + { + new (this) engine(std::forward(args)...); + } + + template + friend bool operator==(const engine&, + const engine&); + + template + friend itype1 operator-(const engine&, + const engine&); + + template + friend std::basic_ostream& + operator<<(std::basic_ostream& out, + const engine&); + + template + friend std::basic_istream& + operator>>(std::basic_istream& in, + engine& rng); +}; + +template +std::basic_ostream& +operator<<(std::basic_ostream& out, + const engine& rng) +{ + auto orig_flags = out.flags(std::ios_base::dec | std::ios_base::left); + auto space = out.widen(' '); + auto orig_fill = out.fill(); + + out << rng.multiplier() << space + << rng.increment() << space + << rng.state_; + + out.flags(orig_flags); + out.fill(orig_fill); + return out; +} + + +template +std::basic_istream& +operator>>(std::basic_istream& in, + engine& rng) +{ + auto orig_flags = in.flags(std::ios_base::dec | std::ios_base::skipws); + + itype multiplier, increment, state; + in >> multiplier >> increment >> state; + + if (!in.fail()) { + bool good = true; + if (multiplier != rng.multiplier()) { + good = false; + } else if (rng.can_specify_stream) { + rng.set_stream(increment >> 1); + } else if (increment != rng.increment()) { + good = false; + } + if (good) { + rng.state_ = state; + } else { + in.clear(std::ios::failbit); + } + } + + in.flags(orig_flags); + return in; +} + + +template +itype engine::advance( + itype state, itype delta, itype cur_mult, itype cur_plus) +{ + // The method used here is based on Brown, "Random Number Generation + // with Arbitrary Stride,", Transactions of the American Nuclear + // Society (Nov. 1994). The algorithm is very similar to fast + // exponentiation. + // + // Even though delta is an unsigned integer, we can pass a + // signed integer to go backwards, it just goes "the long way round". + + constexpr itype ZERO = 0u; // itype may be a non-trivial types, so + constexpr itype ONE = 1u; // we define some ugly constants. + itype acc_mult = 1; + itype acc_plus = 0; + while (delta > ZERO) { + if (delta & ONE) { + acc_mult *= cur_mult; + acc_plus = acc_plus*cur_mult + cur_plus; + } + cur_plus = (cur_mult+ONE)*cur_plus; + cur_mult *= cur_mult; + delta >>= 1; + } + return acc_mult * state + acc_plus; +} + +template +itype engine::distance( + itype cur_state, itype newstate, itype cur_mult, itype cur_plus, itype mask) +{ + constexpr itype ONE = 1u; // itype could be weird, so use constant + itype the_bit = stream_mixin::is_mcg ? itype(4u) : itype(1u); + itype distance = 0u; + while ((cur_state & mask) != (newstate & mask)) { + if ((cur_state & the_bit) != (newstate & the_bit)) { + cur_state = cur_state * cur_mult + cur_plus; + distance |= the_bit; + } + assert((cur_state & the_bit) == (newstate & the_bit)); + the_bit <<= 1; + cur_plus = (cur_mult+ONE)*cur_plus; + cur_mult *= cur_mult; + } + return stream_mixin::is_mcg ? distance >> 2 : distance; +} + +template +itype operator-(const engine& lhs, + const engine& rhs) +{ + if (lhs.multiplier() != rhs.multiplier() + || lhs.increment() != rhs.increment()) + throw std::logic_error("incomparable generators"); + return rhs.distance(lhs.state_); +} + + +template +bool operator==(const engine& lhs, + const engine& rhs) +{ + return (lhs.multiplier() == rhs.multiplier()) + && (lhs.increment() == rhs.increment()) + && (lhs.state_ == rhs.state_); +} + +template +inline bool operator!=(const engine& lhs, + const engine& rhs) +{ + return !operator==(lhs,rhs); +} + + +template class output_mixin, + bool output_previous = (sizeof(itype) <= 8)> +using oneseq_base = engine, output_previous, + oneseq_stream >; + +template class output_mixin, + bool output_previous = (sizeof(itype) <= 8)> +using unique_base = engine, output_previous, + unique_stream >; + +template class output_mixin, + bool output_previous = (sizeof(itype) <= 8)> +using setseq_base = engine, output_previous, + specific_stream >; + +template class output_mixin, + bool output_previous = (sizeof(itype) <= 8)> +using mcg_base = engine, output_previous, + no_stream >; + +/* + * OUTPUT FUNCTIONS. + * + * These are the core of the PCG generation scheme. They specify how to + * turn the base LCG's internal state into the output value of the final + * generator. + * + * They're implemented as mixin classes. + * + * All of the classes have code that is written to allow it to be applied + * at *arbitrary* bit sizes, although in practice they'll only be used at + * standard sizes supported by C++. + */ + +/* + * XSH RS -- high xorshift, followed by a random shift + * + * Fast. A good performer. + */ + +template +struct xsh_rs_mixin { + static xtype output(itype internal) + { + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8); + constexpr bitcount_t sparebits = bits - xtypebits; + constexpr bitcount_t opbits = + sparebits-5 >= 64 ? 5 + : sparebits-4 >= 32 ? 4 + : sparebits-3 >= 16 ? 3 + : sparebits-2 >= 4 ? 2 + : sparebits-1 >= 1 ? 1 + : 0; + constexpr bitcount_t mask = (1 << opbits) - 1; + constexpr bitcount_t maxrandshift = mask; + constexpr bitcount_t topspare = opbits; + constexpr bitcount_t bottomspare = sparebits - topspare; + constexpr bitcount_t xshift = topspare + (xtypebits+maxrandshift)/2; + bitcount_t rshift = + opbits ? bitcount_t(internal >> (bits - opbits)) & mask : 0; + internal ^= internal >> xshift; + xtype result = xtype(internal >> (bottomspare - maxrandshift + rshift)); + return result; + } +}; + +/* + * XSH RR -- high xorshift, followed by a random rotate + * + * Fast. A good performer. Slightly better statistically than XSH RS. + */ + +template +struct xsh_rr_mixin { + static xtype output(itype internal) + { + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype)*8); + constexpr bitcount_t sparebits = bits - xtypebits; + constexpr bitcount_t wantedopbits = + xtypebits >= 128 ? 7 + : xtypebits >= 64 ? 6 + : xtypebits >= 32 ? 5 + : xtypebits >= 16 ? 4 + : 3; + constexpr bitcount_t opbits = + sparebits >= wantedopbits ? wantedopbits + : sparebits; + constexpr bitcount_t amplifier = wantedopbits - opbits; + constexpr bitcount_t mask = (1 << opbits) - 1; + constexpr bitcount_t topspare = opbits; + constexpr bitcount_t bottomspare = sparebits - topspare; + constexpr bitcount_t xshift = (topspare + xtypebits)/2; + bitcount_t rot = opbits ? bitcount_t(internal >> (bits - opbits)) & mask + : 0; + bitcount_t amprot = (rot << amplifier) & mask; + internal ^= internal >> xshift; + xtype result = xtype(internal >> bottomspare); + result = rotr(result, amprot); + return result; + } +}; + +/* + * RXS -- random xorshift + */ + +template +struct rxs_mixin { +static xtype output_rxs(itype internal) + { + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype)*8); + constexpr bitcount_t shift = bits - xtypebits; + constexpr bitcount_t extrashift = (xtypebits - shift)/2; + bitcount_t rshift = shift > 64+8 ? (internal >> (bits - 6)) & 63 + : shift > 32+4 ? (internal >> (bits - 5)) & 31 + : shift > 16+2 ? (internal >> (bits - 4)) & 15 + : shift > 8+1 ? (internal >> (bits - 3)) & 7 + : shift > 4+1 ? (internal >> (bits - 2)) & 3 + : shift > 2+1 ? (internal >> (bits - 1)) & 1 + : 0; + internal ^= internal >> (shift + extrashift - rshift); + xtype result = internal >> rshift; + return result; + } +}; + +/* + * RXS M XS -- random xorshift, mcg multiply, fixed xorshift + * + * The most statistically powerful generator, but all those steps + * make it slower than some of the others. We give it the rottenest jobs. + * + * Because it's usually used in contexts where the state type and the + * result type are the same, it is a permutation and is thus invertable. + * We thus provide a function to invert it. This function is used to + * for the "inside out" generator used by the extended generator. + */ + +/* Defined type-based concepts for the multiplication step. They're actually + * all derived by truncating the 128-bit, which was computed to be a good + * "universal" constant. + */ + +template +struct mcg_multiplier { + // Not defined for an arbitrary type +}; + +template +struct mcg_unmultiplier { + // Not defined for an arbitrary type +}; + +PCG_DEFINE_CONSTANT(uint8_t, mcg, multiplier, 217U) +PCG_DEFINE_CONSTANT(uint8_t, mcg, unmultiplier, 105U) + +PCG_DEFINE_CONSTANT(uint16_t, mcg, multiplier, 62169U) +PCG_DEFINE_CONSTANT(uint16_t, mcg, unmultiplier, 28009U) + +PCG_DEFINE_CONSTANT(uint32_t, mcg, multiplier, 277803737U) +PCG_DEFINE_CONSTANT(uint32_t, mcg, unmultiplier, 2897767785U) + +PCG_DEFINE_CONSTANT(uint64_t, mcg, multiplier, 12605985483714917081ULL) +PCG_DEFINE_CONSTANT(uint64_t, mcg, unmultiplier, 15009553638781119849ULL) + +PCG_DEFINE_CONSTANT(pcg128_t, mcg, multiplier, + PCG_128BIT_CONSTANT(17766728186571221404ULL, 12605985483714917081ULL)) +PCG_DEFINE_CONSTANT(pcg128_t, mcg, unmultiplier, + PCG_128BIT_CONSTANT(14422606686972528997ULL, 15009553638781119849ULL)) + + +template +struct rxs_m_xs_mixin { + static xtype output(itype internal) + { + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8); + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t opbits = xtypebits >= 128 ? 6 + : xtypebits >= 64 ? 5 + : xtypebits >= 32 ? 4 + : xtypebits >= 16 ? 3 + : 2; + constexpr bitcount_t shift = bits - xtypebits; + constexpr bitcount_t mask = (1 << opbits) - 1; + bitcount_t rshift = + opbits ? bitcount_t(internal >> (bits - opbits)) & mask : 0; + internal ^= internal >> (opbits + rshift); + internal *= mcg_multiplier::multiplier(); + xtype result = internal >> shift; + result ^= result >> ((2U*xtypebits+2U)/3U); + return result; + } + + static itype unoutput(itype internal) + { + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t opbits = bits >= 128 ? 6 + : bits >= 64 ? 5 + : bits >= 32 ? 4 + : bits >= 16 ? 3 + : 2; + constexpr bitcount_t mask = (1 << opbits) - 1; + + internal = unxorshift(internal, bits, (2U*bits+2U)/3U); + + internal *= mcg_unmultiplier::unmultiplier(); + + bitcount_t rshift = opbits ? (internal >> (bits - opbits)) & mask : 0; + internal = unxorshift(internal, bits, opbits + rshift); + + return internal; + } +}; + + +/* + * RXS M -- random xorshift, mcg multiply + */ + +template +struct rxs_m_mixin { + static xtype output(itype internal) + { + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8); + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t opbits = xtypebits >= 128 ? 6 + : xtypebits >= 64 ? 5 + : xtypebits >= 32 ? 4 + : xtypebits >= 16 ? 3 + : 2; + constexpr bitcount_t shift = bits - xtypebits; + constexpr bitcount_t mask = (1 << opbits) - 1; + bitcount_t rshift = opbits ? (internal >> (bits - opbits)) & mask : 0; + internal ^= internal >> (opbits + rshift); + internal *= mcg_multiplier::multiplier(); + xtype result = internal >> shift; + return result; + } +}; + +/* + * XSL RR -- fixed xorshift (to low bits), random rotate + * + * Useful for 128-bit types that are split across two CPU registers. + */ + +template +struct xsl_rr_mixin { + static xtype output(itype internal) + { + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8); + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t sparebits = bits - xtypebits; + constexpr bitcount_t wantedopbits = xtypebits >= 128 ? 7 + : xtypebits >= 64 ? 6 + : xtypebits >= 32 ? 5 + : xtypebits >= 16 ? 4 + : 3; + constexpr bitcount_t opbits = sparebits >= wantedopbits ? wantedopbits + : sparebits; + constexpr bitcount_t amplifier = wantedopbits - opbits; + constexpr bitcount_t mask = (1 << opbits) - 1; + constexpr bitcount_t topspare = sparebits; + constexpr bitcount_t bottomspare = sparebits - topspare; + constexpr bitcount_t xshift = (topspare + xtypebits) / 2; + + bitcount_t rot = + opbits ? bitcount_t(internal >> (bits - opbits)) & mask : 0; + bitcount_t amprot = (rot << amplifier) & mask; + internal ^= internal >> xshift; + xtype result = xtype(internal >> bottomspare); + result = rotr(result, amprot); + return result; + } +}; + + +/* + * XSL RR RR -- fixed xorshift (to low bits), random rotate (both parts) + * + * Useful for 128-bit types that are split across two CPU registers. + * If you really want an invertable 128-bit RNG, I guess this is the one. + */ + +template struct halfsize_trait {}; +template <> struct halfsize_trait { typedef uint64_t type; }; +template <> struct halfsize_trait { typedef uint32_t type; }; +template <> struct halfsize_trait { typedef uint16_t type; }; +template <> struct halfsize_trait { typedef uint8_t type; }; + +template +struct xsl_rr_rr_mixin { + typedef typename halfsize_trait::type htype; + + static itype output(itype internal) + { + constexpr bitcount_t htypebits = bitcount_t(sizeof(htype) * 8); + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t sparebits = bits - htypebits; + constexpr bitcount_t wantedopbits = htypebits >= 128 ? 7 + : htypebits >= 64 ? 6 + : htypebits >= 32 ? 5 + : htypebits >= 16 ? 4 + : 3; + constexpr bitcount_t opbits = sparebits >= wantedopbits ? wantedopbits + : sparebits; + constexpr bitcount_t amplifier = wantedopbits - opbits; + constexpr bitcount_t mask = (1 << opbits) - 1; + constexpr bitcount_t topspare = sparebits; + constexpr bitcount_t xshift = (topspare + htypebits) / 2; + + bitcount_t rot = + opbits ? bitcount_t(internal >> (bits - opbits)) & mask : 0; + bitcount_t amprot = (rot << amplifier) & mask; + internal ^= internal >> xshift; + htype lowbits = htype(internal); + lowbits = rotr(lowbits, amprot); + htype highbits = htype(internal >> topspare); + bitcount_t rot2 = lowbits & mask; + bitcount_t amprot2 = (rot2 << amplifier) & mask; + highbits = rotr(highbits, amprot2); + return (itype(highbits) << topspare) ^ itype(lowbits); + } +}; + + +/* + * XSH -- fixed xorshift (to high bits) + * + * You shouldn't use this at 64-bits or less. + */ + +template +struct xsh_mixin { + static xtype output(itype internal) + { + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8); + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t sparebits = bits - xtypebits; + constexpr bitcount_t topspare = 0; + constexpr bitcount_t bottomspare = sparebits - topspare; + constexpr bitcount_t xshift = (topspare + xtypebits) / 2; + + internal ^= internal >> xshift; + xtype result = internal >> bottomspare; + return result; + } +}; + +/* + * XSL -- fixed xorshift (to low bits) + * + * You shouldn't use this at 64-bits or less. + */ + +template +struct xsl_mixin { + inline xtype output(itype internal) + { + constexpr bitcount_t xtypebits = bitcount_t(sizeof(xtype) * 8); + constexpr bitcount_t bits = bitcount_t(sizeof(itype) * 8); + constexpr bitcount_t sparebits = bits - xtypebits; + constexpr bitcount_t topspare = sparebits; + constexpr bitcount_t bottomspare = sparebits - topspare; + constexpr bitcount_t xshift = (topspare + xtypebits) / 2; + + internal ^= internal >> xshift; + xtype result = internal >> bottomspare; + return result; + } +}; + +/* ---- End of Output Functions ---- */ + + +template +struct inside_out : private baseclass { + inside_out() = delete; + + typedef typename baseclass::result_type result_type; + typedef typename baseclass::state_type state_type; + static_assert(sizeof(result_type) == sizeof(state_type), + "Require a RNG whose output function is a permutation"); + + static bool external_step(result_type& randval, size_t i) + { + state_type state = baseclass::unoutput(randval); + state = state * baseclass::multiplier() + baseclass::increment() + + state_type(i*2); + result_type result = baseclass::output(state); + randval = result; + state_type zero = + baseclass::is_mcg ? state & state_type(3U) : state_type(0U); + return result == zero; + } + + static bool external_advance(result_type& randval, size_t i, + result_type delta, bool forwards = true) + { + state_type state = baseclass::unoutput(randval); + state_type mult = baseclass::multiplier(); + state_type inc = baseclass::increment() + state_type(i*2); + state_type zero = + baseclass::is_mcg ? state & state_type(3U) : state_type(0U); + state_type dist_to_zero = baseclass::distance(state, zero, mult, inc); + bool crosses_zero = + forwards ? dist_to_zero <= delta + : (-dist_to_zero) <= delta; + if (!forwards) + delta = -delta; + state = baseclass::advance(state, delta, mult, inc); + randval = baseclass::output(state); + return crosses_zero; + } +}; + + +template +class extended : public baseclass { +public: + typedef typename baseclass::state_type state_type; + typedef typename baseclass::result_type result_type; + typedef inside_out insideout; + +private: + static constexpr bitcount_t rtypebits = sizeof(result_type)*8; + static constexpr bitcount_t stypebits = sizeof(state_type)*8; + + static constexpr bitcount_t tick_limit_pow2 = 64U; + + static constexpr size_t table_size = 1UL << table_pow2; + static constexpr size_t table_shift = stypebits - table_pow2; + static constexpr state_type table_mask = + (state_type(1U) << table_pow2) - state_type(1U); + + static constexpr bool may_tick = + (advance_pow2 < stypebits) && (advance_pow2 < tick_limit_pow2); + static constexpr size_t tick_shift = stypebits - advance_pow2; + static constexpr state_type tick_mask = + may_tick ? state_type( + (uint64_t(1) << (advance_pow2*may_tick)) - 1) + // ^-- stupidity to appease GCC warnings + : ~state_type(0U); + + static constexpr bool may_tock = stypebits < tick_limit_pow2; + + result_type data_[table_size]; + + PCG_NOINLINE void advance_table(); + + PCG_NOINLINE void advance_table(state_type delta, bool isForwards = true); + + result_type& get_extended_value() + { + state_type state = this->state_; + if (kdd && baseclass::is_mcg) { + // The low order bits of an MCG are constant, so drop them. + state >>= 2; + } + size_t index = kdd ? state & table_mask + : state >> table_shift; + + if (may_tick) { + bool tick = kdd ? (state & tick_mask) == state_type(0u) + : (state >> tick_shift) == state_type(0u); + if (tick) + advance_table(); + } + if (may_tock) { + bool tock = state == state_type(0u); + if (tock) + advance_table(); + } + return data_[index]; + } + +public: + static constexpr size_t period_pow2() + { + return baseclass::period_pow2() + table_size*extvalclass::period_pow2(); + } + + __attribute__((always_inline)) result_type operator()() + { + result_type rhs = get_extended_value(); + result_type lhs = this->baseclass::operator()(); + return lhs ^ rhs; + } + + result_type operator()(result_type upper_bound) + { + return bounded_rand(*this, upper_bound); + } + + void set(result_type wanted) + { + result_type& rhs = get_extended_value(); + result_type lhs = this->baseclass::operator()(); + rhs = lhs ^ wanted; + } + + void advance(state_type distance, bool forwards = true); + + void backstep(state_type distance) + { + advance(distance, false); + } + + extended(const result_type* data) + : baseclass() + { + datainit(data); + } + + extended(const result_type* data, state_type seed) + : baseclass(seed) + { + datainit(data); + } + + // This function may or may not exist. It thus has to be a template + // to use SFINAE; users don't have to worry about its template-ness. + + template + extended(const result_type* data, state_type seed, + typename bc::stream_state stream_seed) + : baseclass(seed, stream_seed) + { + datainit(data); + } + + extended() + : baseclass() + { + selfinit(); + } + + extended(state_type seed) + : baseclass(seed) + { + selfinit(); + } + + // This function may or may not exist. It thus has to be a template + // to use SFINAE; users don't have to worry about its template-ness. + + template + extended(state_type seed, typename bc::stream_state stream_seed) + : baseclass(seed, stream_seed) + { + selfinit(); + } + +private: + void selfinit(); + void datainit(const result_type* data); + +public: + + template::value + && !std::is_convertible::value>::type> + extended(SeedSeq&& seedSeq) + : baseclass(seedSeq) + { + generate_to(seedSeq, data_); + } + + template + void seed(Args&&... args) + { + new (this) extended(std::forward(args)...); + } + + template + friend bool operator==(const extended&, + const extended&); + + template + friend std::basic_ostream& + operator<<(std::basic_ostream& out, + const extended&); + + template + friend std::basic_istream& + operator>>(std::basic_istream& in, + extended&); + +}; + + +template +void extended::datainit( + const result_type* data) +{ + for (size_t i = 0; i < table_size; ++i) + data_[i] = data[i]; +} + +template +void extended::selfinit() +{ + // We need to fill the extended table with something, and we have + // very little provided data, so we use the base generator to + // produce values. Although not ideal (use a seed sequence, folks!), + // unexpected correlations are mitigated by + // - using XOR differences rather than the number directly + // - the way the table is accessed, its values *won't* be accessed + // in the same order the were written. + // - any strange correlations would only be apparent if we + // were to backstep the generator so that the base generator + // was generating the same values again + result_type xdiff = baseclass::operator()() - baseclass::operator()(); + for (size_t i = 0; i < table_size; ++i) { + data_[i] = baseclass::operator()() ^ xdiff; + } +} + +template +bool operator==(const extended& lhs, + const extended& rhs) +{ + auto& base_lhs = static_cast(lhs); + auto& base_rhs = static_cast(rhs); + return base_lhs == base_rhs + && !memcmp((void*) lhs.data_, (void*) rhs.data_, sizeof(lhs.data_)); +} + +template +inline bool operator!=(const extended& lhs, + const extended& rhs) +{ + return lhs != rhs; +} + +template +std::basic_ostream& +operator<<(std::basic_ostream& out, + const extended& rng) +{ + auto orig_flags = out.flags(std::ios_base::dec | std::ios_base::left); + auto space = out.widen(' '); + auto orig_fill = out.fill(); + + out << rng.multiplier() << space + << rng.increment() << space + << rng.state_; + + for (const auto& datum : rng.data_) + out << space << datum; + + out.flags(orig_flags); + out.fill(orig_fill); + return out; +} + +template +std::basic_istream& +operator>>(std::basic_istream& in, + extended& rng) +{ + extended new_rng; + auto& base_rng = static_cast(new_rng); + in >> base_rng; + + if (in.fail()) + return in; + + auto orig_flags = in.flags(std::ios_base::dec | std::ios_base::skipws); + + for (auto& datum : new_rng.data_) { + in >> datum; + if (in.fail()) + goto bail; + } + + rng = new_rng; + +bail: + in.flags(orig_flags); + return in; +} + + + +template +void +extended::advance_table() +{ + bool carry = false; + for (size_t i = 0; i < table_size; ++i) { + if (carry) { + carry = insideout::external_step(data_[i],i+1); + } + bool carry2 = insideout::external_step(data_[i],i+1); + carry = carry || carry2; + } +} + +template +void +extended::advance_table( + state_type delta, bool isForwards) +{ + typedef typename baseclass::state_type base_state_t; + typedef typename extvalclass::state_type ext_state_t; + constexpr bitcount_t basebits = sizeof(base_state_t)*8; + constexpr bitcount_t extbits = sizeof(ext_state_t)*8; + static_assert(basebits <= extbits || advance_pow2 > 0, + "Current implementation might overflow its carry"); + + base_state_t carry = 0; + for (size_t i = 0; i < table_size; ++i) { + base_state_t total_delta = carry + delta; + ext_state_t trunc_delta = ext_state_t(total_delta); + if (basebits > extbits) { + carry = total_delta >> extbits; + } else { + carry = 0; + } + carry += + insideout::external_advance(data_[i],i+1, trunc_delta, isForwards); + } +} + +template +void extended::advance( + state_type distance, bool forwards) +{ + static_assert(kdd, + "Efficient advance is too hard for non-kdd extension. " + "For a weak advance, cast to base class"); + state_type zero = + baseclass::is_mcg ? this->state_ & state_type(3U) : state_type(0U); + if (may_tick) { + state_type ticks = distance >> (advance_pow2*may_tick); + // ^-- stupidity to appease GCC + // warnings + state_type adv_mask = + baseclass::is_mcg ? tick_mask << 2 : tick_mask; + state_type next_advance_distance = this->distance(zero, adv_mask); + if (!forwards) + next_advance_distance = (-next_advance_distance) & tick_mask; + if (next_advance_distance < (distance & tick_mask)) { + ++ticks; + } + if (ticks) + advance_table(ticks, forwards); + } + if (forwards) { + if (may_tock && this->distance(zero) <= distance) + advance_table(); + baseclass::advance(distance); + } else { + if (may_tock && -(this->distance(zero)) <= distance) + advance_table(state_type(1U), false); + baseclass::advance(-distance); + } +} + +} // namespace pcg_detail + +namespace pcg_engines { + +using namespace pcg_detail; + +/* Predefined types for XSH RS */ + +typedef oneseq_base oneseq_xsh_rs_16_8; +typedef oneseq_base oneseq_xsh_rs_32_16; +typedef oneseq_base oneseq_xsh_rs_64_32; +typedef oneseq_base oneseq_xsh_rs_128_64; + +typedef unique_base unique_xsh_rs_16_8; +typedef unique_base unique_xsh_rs_32_16; +typedef unique_base unique_xsh_rs_64_32; +typedef unique_base unique_xsh_rs_128_64; + +typedef setseq_base setseq_xsh_rs_16_8; +typedef setseq_base setseq_xsh_rs_32_16; +typedef setseq_base setseq_xsh_rs_64_32; +typedef setseq_base setseq_xsh_rs_128_64; + +typedef mcg_base mcg_xsh_rs_16_8; +typedef mcg_base mcg_xsh_rs_32_16; +typedef mcg_base mcg_xsh_rs_64_32; +typedef mcg_base mcg_xsh_rs_128_64; + +/* Predefined types for XSH RR */ + +typedef oneseq_base oneseq_xsh_rr_16_8; +typedef oneseq_base oneseq_xsh_rr_32_16; +typedef oneseq_base oneseq_xsh_rr_64_32; +typedef oneseq_base oneseq_xsh_rr_128_64; + +typedef unique_base unique_xsh_rr_16_8; +typedef unique_base unique_xsh_rr_32_16; +typedef unique_base unique_xsh_rr_64_32; +typedef unique_base unique_xsh_rr_128_64; + +typedef setseq_base setseq_xsh_rr_16_8; +typedef setseq_base setseq_xsh_rr_32_16; +typedef setseq_base setseq_xsh_rr_64_32; +typedef setseq_base setseq_xsh_rr_128_64; + +typedef mcg_base mcg_xsh_rr_16_8; +typedef mcg_base mcg_xsh_rr_32_16; +typedef mcg_base mcg_xsh_rr_64_32; +typedef mcg_base mcg_xsh_rr_128_64; + + +/* Predefined types for RXS M XS */ + +typedef oneseq_base oneseq_rxs_m_xs_8_8; +typedef oneseq_base oneseq_rxs_m_xs_16_16; +typedef oneseq_base oneseq_rxs_m_xs_32_32; +typedef oneseq_base oneseq_rxs_m_xs_64_64; +typedef oneseq_base oneseq_rxs_m_xs_128_128; + +typedef unique_base unique_rxs_m_xs_8_8; +typedef unique_base unique_rxs_m_xs_16_16; +typedef unique_base unique_rxs_m_xs_32_32; +typedef unique_base unique_rxs_m_xs_64_64; +typedef unique_base unique_rxs_m_xs_128_128; + +typedef setseq_base setseq_rxs_m_xs_8_8; +typedef setseq_base setseq_rxs_m_xs_16_16; +typedef setseq_base setseq_rxs_m_xs_32_32; +typedef setseq_base setseq_rxs_m_xs_64_64; +typedef setseq_base setseq_rxs_m_xs_128_128; + + // MCG versions don't make sense here, so aren't defined. + +/* Predefined types for XSL RR (only defined for "large" types) */ + +typedef oneseq_base oneseq_xsl_rr_64_32; +typedef oneseq_base oneseq_xsl_rr_128_64; + +typedef unique_base unique_xsl_rr_64_32; +typedef unique_base unique_xsl_rr_128_64; + +typedef setseq_base setseq_xsl_rr_64_32; +typedef setseq_base setseq_xsl_rr_128_64; + +typedef mcg_base mcg_xsl_rr_64_32; +typedef mcg_base mcg_xsl_rr_128_64; + + +/* Predefined types for XSL RR RR (only defined for "large" types) */ + +typedef oneseq_base + oneseq_xsl_rr_rr_64_64; +typedef oneseq_base + oneseq_xsl_rr_rr_128_128; + +typedef unique_base + unique_xsl_rr_rr_64_64; +typedef unique_base + unique_xsl_rr_rr_128_128; + +typedef setseq_base + setseq_xsl_rr_rr_64_64; +typedef setseq_base + setseq_xsl_rr_rr_128_128; + + // MCG versions don't make sense here, so aren't defined. + +/* Extended generators */ + +template +using ext_std8 = extended; + +template +using ext_std16 = extended; + +template +using ext_std32 = extended; + +template +using ext_std64 = extended; + + +template +using ext_oneseq_rxs_m_xs_32_32 = + ext_std32; + +template +using ext_mcg_xsh_rs_64_32 = + ext_std32; + +template +using ext_oneseq_xsh_rs_64_32 = + ext_std32; + +template +using ext_setseq_xsh_rr_64_32 = + ext_std32; + +template +using ext_mcg_xsl_rr_128_64 = + ext_std64; + +template +using ext_oneseq_xsl_rr_128_64 = + ext_std64; + +template +using ext_setseq_xsl_rr_128_64 = + ext_std64; + +} // namespace pcg_engines + +typedef pcg_engines::setseq_xsh_rr_64_32 pcg32; +typedef pcg_engines::oneseq_xsh_rr_64_32 pcg32_oneseq; +typedef pcg_engines::unique_xsh_rr_64_32 pcg32_unique; +typedef pcg_engines::mcg_xsh_rs_64_32 pcg32_fast; + +typedef pcg_engines::setseq_xsl_rr_128_64 pcg64; +typedef pcg_engines::oneseq_xsl_rr_128_64 pcg64_oneseq; +typedef pcg_engines::unique_xsl_rr_128_64 pcg64_unique; +typedef pcg_engines::mcg_xsl_rr_128_64 pcg64_fast; + +typedef pcg_engines::setseq_rxs_m_xs_8_8 pcg8_once_insecure; +typedef pcg_engines::setseq_rxs_m_xs_16_16 pcg16_once_insecure; +typedef pcg_engines::setseq_rxs_m_xs_32_32 pcg32_once_insecure; +typedef pcg_engines::setseq_rxs_m_xs_64_64 pcg64_once_insecure; +typedef pcg_engines::setseq_xsl_rr_rr_128_128 pcg128_once_insecure; + +typedef pcg_engines::oneseq_rxs_m_xs_8_8 pcg8_oneseq_once_insecure; +typedef pcg_engines::oneseq_rxs_m_xs_16_16 pcg16_oneseq_once_insecure; +typedef pcg_engines::oneseq_rxs_m_xs_32_32 pcg32_oneseq_once_insecure; +typedef pcg_engines::oneseq_rxs_m_xs_64_64 pcg64_oneseq_once_insecure; +typedef pcg_engines::oneseq_xsl_rr_rr_128_128 pcg128_oneseq_once_insecure; + + +// These two extended RNGs provide two-dimensionally equidistributed +// 32-bit generators. pcg32_k2_fast occupies the same space as pcg64, +// and can be called twice to generate 64 bits, but does not required +// 128-bit math; on 32-bit systems, it's faster than pcg64 as well. + +typedef pcg_engines::ext_setseq_xsh_rr_64_32<6,16,true> pcg32_k2; +typedef pcg_engines::ext_oneseq_xsh_rs_64_32<6,32,true> pcg32_k2_fast; + +// These eight extended RNGs have about as much state as arc4random +// +// - the k variants are k-dimensionally equidistributed +// - the c variants offer better crypographic security +// +// (just how good the cryptographic security is is an open question) + +typedef pcg_engines::ext_setseq_xsh_rr_64_32<6,16,true> pcg32_k64; +typedef pcg_engines::ext_mcg_xsh_rs_64_32<6,32,true> pcg32_k64_oneseq; +typedef pcg_engines::ext_oneseq_xsh_rs_64_32<6,32,true> pcg32_k64_fast; + +typedef pcg_engines::ext_setseq_xsh_rr_64_32<6,16,false> pcg32_c64; +typedef pcg_engines::ext_oneseq_xsh_rs_64_32<6,32,false> pcg32_c64_oneseq; +typedef pcg_engines::ext_mcg_xsh_rs_64_32<6,32,false> pcg32_c64_fast; + +typedef pcg_engines::ext_setseq_xsl_rr_128_64<5,16,true> pcg64_k32; +typedef pcg_engines::ext_oneseq_xsl_rr_128_64<5,128,true> pcg64_k32_oneseq; +typedef pcg_engines::ext_mcg_xsl_rr_128_64<5,128,true> pcg64_k32_fast; + +typedef pcg_engines::ext_setseq_xsl_rr_128_64<5,16,false> pcg64_c32; +typedef pcg_engines::ext_oneseq_xsl_rr_128_64<5,128,false> pcg64_c32_oneseq; +typedef pcg_engines::ext_mcg_xsl_rr_128_64<5,128,false> pcg64_c32_fast; + +// These eight extended RNGs have more state than the Mersenne twister +// +// - the k variants are k-dimensionally equidistributed +// - the c variants offer better crypographic security +// +// (just how good the cryptographic security is is an open question) + +typedef pcg_engines::ext_setseq_xsh_rr_64_32<10,16,true> pcg32_k1024; +typedef pcg_engines::ext_oneseq_xsh_rs_64_32<10,32,true> pcg32_k1024_fast; + +typedef pcg_engines::ext_setseq_xsh_rr_64_32<10,16,false> pcg32_c1024; +typedef pcg_engines::ext_oneseq_xsh_rs_64_32<10,32,false> pcg32_c1024_fast; + +typedef pcg_engines::ext_setseq_xsl_rr_128_64<10,16,true> pcg64_k1024; +typedef pcg_engines::ext_oneseq_xsl_rr_128_64<10,128,true> pcg64_k1024_fast; + +typedef pcg_engines::ext_setseq_xsl_rr_128_64<10,16,false> pcg64_c1024; +typedef pcg_engines::ext_oneseq_xsl_rr_128_64<10,128,false> pcg64_c1024_fast; + +// These generators have an insanely huge period (2^524352), and is suitable +// for silly party tricks, such as dumping out 64 KB ZIP files at an arbitrary +// point in the future. [Actually, over the full period of the generator, it +// will produce every 64 KB ZIP file 2^64 times!] + +typedef pcg_engines::ext_setseq_xsh_rr_64_32<14,16,true> pcg32_k16384; +typedef pcg_engines::ext_oneseq_xsh_rs_64_32<14,32,true> pcg32_k16384_fast; + +#endif // PCG_RAND_HPP_INCLUDED diff --git a/src/pcg_uint128.hpp b/src/pcg_uint128.hpp new file mode 100644 index 0000000000..a0ad954ef5 --- /dev/null +++ b/src/pcg_uint128.hpp @@ -0,0 +1,750 @@ +/* + * PCG Random Number Generation for C++ + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code provides a a C++ class that can provide 128-bit (or higher) + * integers. To produce 2K-bit integers, it uses two K-bit integers, + * placed in a union that allowes the code to also see them as four K/2 bit + * integers (and access them either directly name, or by index). + * + * It may seem like we're reinventing the wheel here, because several + * libraries already exist that support large integers, but most existing + * libraries provide a very generic multiprecision code, but here we're + * operating at a fixed size. Also, most other libraries are fairly + * heavyweight. So we use a direct implementation. Sadly, it's much slower + * than hand-coded assembly or direct CPU support. + */ + +#ifndef PCG_UINT128_HPP_INCLUDED +#define PCG_UINT128_HPP_INCLUDED 1 + +#include +#include +#include +#include +#include +#include +#include + +/* + * We want to lay the type out the same way that a native type would be laid + * out, which means we must know the machine's endian, at compile time. + * This ugliness attempts to do so. + */ + +#ifndef PCG_LITTLE_ENDIAN + #if defined(__BYTE_ORDER__) + #if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ + #define PCG_LITTLE_ENDIAN 1 + #elif __BYTE_ORDER__ == __ORDER_BIG_ENDIAN__ + #define PCG_LITTLE_ENDIAN 0 + #else + #error __BYTE_ORDER__ does not match a standard endian, pick a side + #endif + #elif __LITTLE_ENDIAN__ || _LITTLE_ENDIAN + #define PCG_LITTLE_ENDIAN 1 + #elif __BIG_ENDIAN__ || _BIG_ENDIAN + #define PCG_LITTLE_ENDIAN 0 + #elif __x86_64 || __x86_64__ || __i386 || __i386__ || _M_IX86 || _M_X64 + #define PCG_LITTLE_ENDIAN 1 + #elif __powerpc__ || __POWERPC__ || __ppc__ || __PPC__ \ + || __m68k__ || __mc68000__ + #define PCG_LITTLE_ENDIAN 0 + #else + #error Unable to determine target endianness + #endif +#endif + +namespace pcg_extras { + +// Recent versions of GCC have intrinsics we can use to quickly calculate +// the number of leading and trailing zeros in a number. If possible, we +// use them, otherwise we fall back to old-fashioned bit twiddling to figure +// them out. + +#ifndef PCG_BITCOUNT_T + typedef uint8_t bitcount_t; +#else + typedef PCG_BITCOUNT_T bitcount_t; +#endif + +/* + * Provide some useful helper functions + * * flog2 floor(log2(x)) + * * trailingzeros number of trailing zero bits + */ + +#ifdef __GNUC__ // Any GNU-compatible compiler supporting C++11 has + // some useful intrinsics we can use. + +inline bitcount_t flog2(uint32_t v) +{ + return 31 - __builtin_clz(v); +} + +inline bitcount_t trailingzeros(uint32_t v) +{ + return __builtin_ctz(v); +} + +inline bitcount_t flog2(uint64_t v) +{ +#if UINT64_MAX == ULONG_MAX + return 63 - __builtin_clzl(v); +#elif UINT64_MAX == ULLONG_MAX + return 63 - __builtin_clzll(v); +#else + #error Cannot find a function for uint64_t +#endif +} + +inline bitcount_t trailingzeros(uint64_t v) +{ +#if UINT64_MAX == ULONG_MAX + return __builtin_ctzl(v); +#elif UINT64_MAX == ULLONG_MAX + return __builtin_ctzll(v); +#else + #error Cannot find a function for uint64_t +#endif +} + +#else // Otherwise, we fall back to bit twiddling + // implementations + +inline bitcount_t flog2(uint32_t v) +{ + // Based on code by Eric Cole and Mark Dickinson, which appears at + // https://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn + + static const uint8_t multiplyDeBruijnBitPos[32] = { + 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, + 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31 + }; + + v |= v >> 1; // first round down to one less than a power of 2 + v |= v >> 2; + v |= v >> 4; + v |= v >> 8; + v |= v >> 16; + + return multiplyDeBruijnBitPos[(uint32_t)(v * 0x07C4ACDDU) >> 27]; +} + +inline bitcount_t trailingzeros(uint32_t v) +{ + static const uint8_t multiplyDeBruijnBitPos[32] = { + 0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, + 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9 + }; + + return multiplyDeBruijnBitPos[((uint32_t)((v & -v) * 0x077CB531U)) >> 27]; +} + +inline bitcount_t flog2(uint64_t v) +{ + uint32_t high = v >> 32; + uint32_t low = uint32_t(v); + + return high ? 32+flog2(high) : flog2(low); +} + +inline bitcount_t trailingzeros(uint64_t v) +{ + uint32_t high = v >> 32; + uint32_t low = uint32_t(v); + + return low ? trailingzeros(low) : trailingzeros(high)+32; +} + +#endif + +template +inline bitcount_t clog2(UInt v) +{ + return flog2(v) + ((v & (-v)) != v); +} + +template +inline UInt addwithcarry(UInt x, UInt y, bool carryin, bool* carryout) +{ + UInt half_result = y + carryin; + UInt result = x + half_result; + *carryout = (half_result < y) || (result < x); + return result; +} + +template +inline UInt subwithcarry(UInt x, UInt y, bool carryin, bool* carryout) +{ + UInt half_result = y + carryin; + UInt result = x - half_result; + *carryout = (half_result < y) || (result > x); + return result; +} + + +template +class uint_x4 { +// private: +public: + union { +#if PCG_LITTLE_ENDIAN + struct { + UInt v0, v1, v2, v3; + } w; + struct { + UIntX2 v01, v23; + } d; +#else + struct { + UInt v3, v2, v1, v0; + } w; + struct { + UIntX2 v23, v01; + } d; +#endif + // For the array access versions, the code that uses the array + // must handle endian itself. Yuck. + UInt wa[4]; + UIntX2 da[2]; + }; + +public: + uint_x4() = default; + + constexpr uint_x4(UInt v3, UInt v2, UInt v1, UInt v0) +#if PCG_LITTLE_ENDIAN + : w{v0, v1, v2, v3} +#else + : w{v3, v2, v1, v0} +#endif + { + // Nothing (else) to do + } + + constexpr uint_x4(UIntX2 v23, UIntX2 v01) +#if PCG_LITTLE_ENDIAN + : d{v01,v23} +#else + : d{v23,v01} +#endif + { + // Nothing (else) to do + } + + template::value + && sizeof(Integral) <= sizeof(UIntX2)) + >::type* = nullptr> + constexpr uint_x4(Integral v01) +#if PCG_LITTLE_ENDIAN + : d{UIntX2(v01),0UL} +#else + : d{0UL,UIntX2(v01)} +#endif + { + // Nothing (else) to do + } + + explicit constexpr operator uint64_t() const + { + return d.v01; + } + + explicit constexpr operator uint32_t() const + { + return w.v0; + } + + explicit constexpr operator int() const + { + return w.v0; + } + + explicit constexpr operator uint16_t() const + { + return w.v0; + } + + explicit constexpr operator uint8_t() const + { + return w.v0; + } + + typedef typename std::conditional::value, + unsigned long long, + unsigned long>::type + uint_missing_t; + + explicit constexpr operator uint_missing_t() const + { + return d.v01; + } + + explicit constexpr operator bool() const + { + return d.v01 || d.v23; + } + + template + friend uint_x4 operator*(const uint_x4&, const uint_x4&); + + template + friend std::pair< uint_x4,uint_x4 > + divmod(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator+(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator-(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator<<(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator>>(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator&(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator|(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator^(const uint_x4&, const uint_x4&); + + template + friend bool operator==(const uint_x4&, const uint_x4&); + + template + friend bool operator!=(const uint_x4&, const uint_x4&); + + template + friend bool operator<(const uint_x4&, const uint_x4&); + + template + friend bool operator<=(const uint_x4&, const uint_x4&); + + template + friend bool operator>(const uint_x4&, const uint_x4&); + + template + friend bool operator>=(const uint_x4&, const uint_x4&); + + template + friend uint_x4 operator~(const uint_x4&); + + template + friend uint_x4 operator-(const uint_x4&); + + template + friend bitcount_t flog2(const uint_x4&); + + template + friend bitcount_t trailingzeros(const uint_x4&); + + uint_x4& operator*=(const uint_x4& rhs) + { + uint_x4 result = *this * rhs; + return *this = result; + } + + uint_x4& operator/=(const uint_x4& rhs) + { + uint_x4 result = *this / rhs; + return *this = result; + } + + uint_x4& operator%=(const uint_x4& rhs) + { + uint_x4 result = *this % rhs; + return *this = result; + } + + uint_x4& operator+=(const uint_x4& rhs) + { + uint_x4 result = *this + rhs; + return *this = result; + } + + uint_x4& operator-=(const uint_x4& rhs) + { + uint_x4 result = *this - rhs; + return *this = result; + } + + uint_x4& operator&=(const uint_x4& rhs) + { + uint_x4 result = *this & rhs; + return *this = result; + } + + uint_x4& operator|=(const uint_x4& rhs) + { + uint_x4 result = *this | rhs; + return *this = result; + } + + uint_x4& operator^=(const uint_x4& rhs) + { + uint_x4 result = *this ^ rhs; + return *this = result; + } + + uint_x4& operator>>=(bitcount_t shift) + { + uint_x4 result = *this >> shift; + return *this = result; + } + + uint_x4& operator<<=(bitcount_t shift) + { + uint_x4 result = *this << shift; + return *this = result; + } + +}; + +template +bitcount_t flog2(const uint_x4& v) +{ +#if PCG_LITTLE_ENDIAN + for (uint8_t i = 4; i !=0; /* dec in loop */) { + --i; +#else + for (uint8_t i = 0; i < 4; ++i) { +#endif + if (v.wa[i] == 0) + continue; + return flog2(v.wa[i]) + (sizeof(U)*CHAR_BIT)*i; + } + abort(); +} + +template +bitcount_t trailingzeros(const uint_x4& v) +{ +#if PCG_LITTLE_ENDIAN + for (uint8_t i = 0; i < 4; ++i) { +#else + for (uint8_t i = 4; i !=0; /* dec in loop */) { + --i; +#endif + if (v.wa[i] != 0) + return trailingzeros(v.wa[i]) + (sizeof(U)*CHAR_BIT)*i; + } + return (sizeof(U)*CHAR_BIT)*4; +} + +template +std::pair< uint_x4, uint_x4 > + divmod(const uint_x4& orig_dividend, + const uint_x4& divisor) +{ + // If the dividend is less than the divisor, the answer is always zero. + // This takes care of boundary cases like 0/x (which would otherwise be + // problematic because we can't take the log of zero. (The boundary case + // of division by zero is undefined.) + if (orig_dividend < divisor) + return { uint_x4(0UL), orig_dividend }; + + auto dividend = orig_dividend; + + auto log2_divisor = flog2(divisor); + auto log2_dividend = flog2(dividend); + // assert(log2_dividend >= log2_divisor); + bitcount_t logdiff = log2_dividend - log2_divisor; + + constexpr uint_x4 ONE(1UL); + if (logdiff == 0) + return { ONE, dividend - divisor }; + + // Now we change the log difference to + // floor(log2(divisor)) - ceil(log2(dividend)) + // to ensure that we *underestimate* the result. + logdiff -= 1; + + uint_x4 quotient(0UL); + + auto qfactor = ONE << logdiff; + auto factor = divisor << logdiff; + + do { + dividend -= factor; + quotient += qfactor; + while (dividend < factor) { + factor >>= 1; + qfactor >>= 1; + } + } while (dividend >= divisor); + + return { quotient, dividend }; +} + +template +uint_x4 operator/(const uint_x4& dividend, + const uint_x4& divisor) +{ + return divmod(dividend, divisor).first; +} + +template +uint_x4 operator%(const uint_x4& dividend, + const uint_x4& divisor) +{ + return divmod(dividend, divisor).second; +} + + +template +uint_x4 operator*(const uint_x4& a, + const uint_x4& b) +{ + uint_x4 r = {0U, 0U, 0U, 0U}; + bool carryin = false; + bool carryout; + UIntX2 a0b0 = UIntX2(a.w.v0) * UIntX2(b.w.v0); + r.w.v0 = UInt(a0b0); + r.w.v1 = UInt(a0b0 >> 32); + + UIntX2 a1b0 = UIntX2(a.w.v1) * UIntX2(b.w.v0); + r.w.v2 = UInt(a1b0 >> 32); + r.w.v1 = addwithcarry(r.w.v1, UInt(a1b0), carryin, &carryout); + carryin = carryout; + r.w.v2 = addwithcarry(r.w.v2, UInt(0U), carryin, &carryout); + carryin = carryout; + r.w.v3 = addwithcarry(r.w.v3, UInt(0U), carryin, &carryout); + + UIntX2 a0b1 = UIntX2(a.w.v0) * UIntX2(b.w.v1); + carryin = false; + r.w.v2 = addwithcarry(r.w.v2, UInt(a0b1 >> 32), carryin, &carryout); + carryin = carryout; + r.w.v3 = addwithcarry(r.w.v3, UInt(0U), carryin, &carryout); + + carryin = false; + r.w.v1 = addwithcarry(r.w.v1, UInt(a0b1), carryin, &carryout); + carryin = carryout; + r.w.v2 = addwithcarry(r.w.v2, UInt(0U), carryin, &carryout); + carryin = carryout; + r.w.v3 = addwithcarry(r.w.v3, UInt(0U), carryin, &carryout); + + UIntX2 a1b1 = UIntX2(a.w.v1) * UIntX2(b.w.v1); + carryin = false; + r.w.v2 = addwithcarry(r.w.v2, UInt(a1b1), carryin, &carryout); + carryin = carryout; + r.w.v3 = addwithcarry(r.w.v3, UInt(a1b1 >> 32), carryin, &carryout); + + r.d.v23 += a.d.v01 * b.d.v23 + a.d.v23 * b.d.v01; + + return r; +} + + +template +uint_x4 operator+(const uint_x4& a, + const uint_x4& b) +{ + uint_x4 r = {0U, 0U, 0U, 0U}; + + bool carryin = false; + bool carryout; + r.w.v0 = addwithcarry(a.w.v0, b.w.v0, carryin, &carryout); + carryin = carryout; + r.w.v1 = addwithcarry(a.w.v1, b.w.v1, carryin, &carryout); + carryin = carryout; + r.w.v2 = addwithcarry(a.w.v2, b.w.v2, carryin, &carryout); + carryin = carryout; + r.w.v3 = addwithcarry(a.w.v3, b.w.v3, carryin, &carryout); + + return r; +} + +template +uint_x4 operator-(const uint_x4& a, + const uint_x4& b) +{ + uint_x4 r = {0U, 0U, 0U, 0U}; + + bool carryin = false; + bool carryout; + r.w.v0 = subwithcarry(a.w.v0, b.w.v0, carryin, &carryout); + carryin = carryout; + r.w.v1 = subwithcarry(a.w.v1, b.w.v1, carryin, &carryout); + carryin = carryout; + r.w.v2 = subwithcarry(a.w.v2, b.w.v2, carryin, &carryout); + carryin = carryout; + r.w.v3 = subwithcarry(a.w.v3, b.w.v3, carryin, &carryout); + + return r; +} + + +template +uint_x4 operator&(const uint_x4& a, + const uint_x4& b) +{ + return uint_x4(a.d.v23 & b.d.v23, a.d.v01 & b.d.v01); +} + +template +uint_x4 operator|(const uint_x4& a, + const uint_x4& b) +{ + return uint_x4(a.d.v23 | b.d.v23, a.d.v01 | b.d.v01); +} + +template +uint_x4 operator^(const uint_x4& a, + const uint_x4& b) +{ + return uint_x4(a.d.v23 ^ b.d.v23, a.d.v01 ^ b.d.v01); +} + +template +uint_x4 operator~(const uint_x4& v) +{ + return uint_x4(~v.d.v23, ~v.d.v01); +} + +template +uint_x4 operator-(const uint_x4& v) +{ + return uint_x4(0UL,0UL) - v; +} + +template +bool operator==(const uint_x4& a, const uint_x4& b) +{ + return (a.d.v01 == b.d.v01) && (a.d.v23 == b.d.v23); +} + +template +bool operator!=(const uint_x4& a, const uint_x4& b) +{ + return !operator==(a,b); +} + + +template +bool operator<(const uint_x4& a, const uint_x4& b) +{ + return (a.d.v23 < b.d.v23) + || ((a.d.v23 == b.d.v23) && (a.d.v01 < b.d.v01)); +} + +template +bool operator>(const uint_x4& a, const uint_x4& b) +{ + return operator<(b,a); +} + +template +bool operator<=(const uint_x4& a, const uint_x4& b) +{ + return !(operator<(b,a)); +} + +template +bool operator>=(const uint_x4& a, const uint_x4& b) +{ + return !(operator<(a,b)); +} + + + +template +uint_x4 operator<<(const uint_x4& v, + const bitcount_t shift) +{ + uint_x4 r = {0U, 0U, 0U, 0U}; + const bitcount_t bits = sizeof(UInt) * CHAR_BIT; + const bitcount_t bitmask = bits - 1; + const bitcount_t shiftdiv = shift / bits; + const bitcount_t shiftmod = shift & bitmask; + + if (shiftmod) { + UInt carryover = 0; +#if PCG_LITTLE_ENDIAN + for (uint8_t out = shiftdiv, in = 0; out < 4; ++out, ++in) { +#else + for (uint8_t out = 4-shiftdiv, in = 4; out != 0; /* dec in loop */) { + --out, --in; +#endif + r.wa[out] = (v.wa[in] << shiftmod) | carryover; + carryover = (v.wa[in] >> (bits - shiftmod)); + } + } else { +#if PCG_LITTLE_ENDIAN + for (uint8_t out = shiftdiv, in = 0; out < 4; ++out, ++in) { +#else + for (uint8_t out = 4-shiftdiv, in = 4; out != 0; /* dec in loop */) { + --out, --in; +#endif + r.wa[out] = v.wa[in]; + } + } + + return r; +} + +template +uint_x4 operator>>(const uint_x4& v, + const bitcount_t shift) +{ + uint_x4 r = {0U, 0U, 0U, 0U}; + const bitcount_t bits = sizeof(UInt) * CHAR_BIT; + const bitcount_t bitmask = bits - 1; + const bitcount_t shiftdiv = shift / bits; + const bitcount_t shiftmod = shift & bitmask; + + if (shiftmod) { + UInt carryover = 0; +#if PCG_LITTLE_ENDIAN + for (uint8_t out = 4-shiftdiv, in = 4; out != 0; /* dec in loop */) { + --out, --in; +#else + for (uint8_t out = shiftdiv, in = 0; out < 4; ++out, ++in) { +#endif + r.wa[out] = (v.wa[in] >> shiftmod) | carryover; + carryover = (v.wa[in] << (bits - shiftmod)); + } + } else { +#if PCG_LITTLE_ENDIAN + for (uint8_t out = 4-shiftdiv, in = 4; out != 0; /* dec in loop */) { + --out, --in; +#else + for (uint8_t out = shiftdiv, in = 0; out < 4; ++out, ++in) { +#endif + r.wa[out] = v.wa[in]; + } + } + + return r; +} + +} // namespace pcg_extras + +#endif // PCG_UINT128_HPP_INCLUDED diff --git a/src/xoshiro128plusplus.cpp b/src/xoshiro128plusplus.cpp new file mode 100644 index 0000000000..cc24d6a3c9 --- /dev/null +++ b/src/xoshiro128plusplus.cpp @@ -0,0 +1,118 @@ +// NOTE: Original URL: http://prng.di.unimi.it/xoshiro128plusplus.c +/** The Interface code (at the very bottom) was written by DRR. */ +/* Written in 2019 by David Blackman and Sebastiano Vigna (vigna@acm.org) + +To the extent possible under law, the author has dedicated all copyright +and related and neighboring rights to this software to the public domain +worldwide. This software is distributed without any warranty. + +See . */ + +#include + +/* This is xoshiro128++ 1.0, one of our 32-bit all-purpose, rock-solid + generators. It has excellent speed, a state size (128 bits) that is + large enough for mild parallelism, and it passes all tests we are aware + of. + + For generating just single-precision (i.e., 32-bit) floating-point + numbers, xoshiro128+ is even faster. + + The state must be seeded so that it is not everywhere zero. */ + + +static inline uint32_t rotl(const uint32_t x, int k) { + return (x << k) | (x >> (32 - k)); +} + + +static uint32_t s[4]; + +uint32_t next(void) { + const uint32_t result = rotl(s[0] + s[3], 7) + s[0]; + + const uint32_t t = s[1] << 9; + + s[2] ^= s[0]; + s[3] ^= s[1]; + s[1] ^= s[2]; + s[0] ^= s[3]; + + s[2] ^= t; + + s[3] = rotl(s[3], 11); + + return result; +} + + +/* This is the jump function for the generator. It is equivalent + to 2^64 calls to next(); it can be used to generate 2^64 + non-overlapping subsequences for parallel computations. */ + +void jump(void) { + static const uint32_t JUMP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b }; + + uint32_t s0 = 0; + uint32_t s1 = 0; + uint32_t s2 = 0; + uint32_t s3 = 0; + for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) + for(int b = 0; b < 32; b++) { + if (JUMP[i] & UINT32_C(1) << b) { + s0 ^= s[0]; + s1 ^= s[1]; + s2 ^= s[2]; + s3 ^= s[3]; + } + next(); + } + + s[0] = s0; + s[1] = s1; + s[2] = s2; + s[3] = s3; +} + + +/* This is the long-jump function for the generator. It is equivalent to + 2^96 calls to next(); it can be used to generate 2^32 starting points, + from each of which jump() will generate 2^32 non-overlapping + subsequences for parallel distributed computations. */ + +void long_jump(void) { + static const uint32_t LONG_JUMP[] = { 0xb523952e, 0x0b6f099f, 0xccf5a0ef, 0x1c580662 }; + + uint32_t s0 = 0; + uint32_t s1 = 0; + uint32_t s2 = 0; + uint32_t s3 = 0; + for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++) + for(int b = 0; b < 32; b++) { + if (LONG_JUMP[i] & UINT32_C(1) << b) { + s0 ^= s[0]; + s1 ^= s[1]; + s2 ^= s[2]; + s3 ^= s[3]; + } + next(); + } + + s[0] = s0; + s[1] = s1; + s[2] = s2; + s[3] = s3; +} + +// ----- Interface ------------------------------- +void x128pp_seed(int seedIn) { + s[0] = seedIn; + s[1] = seedIn+1; + s[2] = seedIn+2; + s[3] = seedIn+3; +} + +unsigned int x128pp_next() { + return (unsigned int)next(); +} + diff --git a/src/xoshiro128plusplus.h b/src/xoshiro128plusplus.h new file mode 100644 index 0000000000..0f87c07e33 --- /dev/null +++ b/src/xoshiro128plusplus.h @@ -0,0 +1,8 @@ +#ifndef INC_XOSHIRO128PLUSPLUS_H +#define INC_XOSHIRO128PLUSPLUS_H + +void x128pp_seed(int); + +unsigned int x128pp_next(); + +#endif diff --git a/test/Makefile b/test/Makefile index be14162b5c..8fa27ac28f 100644 --- a/test/Makefile +++ b/test/Makefile @@ -477,6 +477,9 @@ test.graft: test.flatten: @-cd Test_Flatten && ./RunTest.sh $(OPT) +test.random: + @-cd Test_Random && ./RunTest.sh $(OPT) + # Every test target should go here. COMPLETETESTS=test.general \ test.strip \ @@ -629,7 +632,8 @@ COMPLETETESTS=test.general \ test.emin \ test.evalplateau \ test.graft \ - test.flatten + test.flatten \ + test.random test.all: $(MAKE) test.complete summary diff --git a/test/MasterTest.sh b/test/MasterTest.sh index 6077ee7098..669c8579ce 100644 --- a/test/MasterTest.sh +++ b/test/MasterTest.sh @@ -32,10 +32,12 @@ # DIFFOPTS : Additional options to pass to CPPTRAJ_DIFF # CPPTRAJ_PROFILE : If 1, end of test profiling with gprof performed. # CPPTRAJ_LONG_TEST : If 1, enable long tests +# CPPTRAJ_RNG : Hold the --rng arg for CPPTRAJ. # Cpptraj binary characteristics # CPPTRAJ_ZLIB : If set CPPTRAJ has zlib support. # CPPTRAJ_BZLIB : If set CPPTRAJ has bzip support. # CPPTRAJ_NETCDFLIB : If set CPPTRAJ has NetCDF support. +# CPPTRAJ_C11_SUPPORT : If set CPPTRAJ has C++11 support. # CPPTRAJ_LIBPME : If set CPPTRAJ was compiled with libPME. # CPPTRAJ_MPILIB : If set CPPTRAJ has MPI support. # CPPTRAJ_MATHLIB : If set CPPTRAJ was compiled with math libraries. @@ -501,7 +503,7 @@ RunCpptraj() { if [ -z "$CPPTRAJ_DACDIF" ] ; then OUT " CPPTRAJ: $1" fi - cpptraj_cmd="$CPPTRAJ_TIME $DO_PARALLEL $VALGRIND $CPPTRAJ $TOP $INPUT $CPPTRAJ_DEBUG" + cpptraj_cmd="$CPPTRAJ_TIME $DO_PARALLEL $VALGRIND $CPPTRAJ $TOP $INPUT $CPPTRAJ_RNG $CPPTRAJ_DEBUG" if [ ! -z "$CPPTRAJ_DEBUG" ] ; then echo "$cpptraj_cmd >> $CPPTRAJ_OUTPUT 2>>$CPPTRAJ_ERROR" fi @@ -750,6 +752,7 @@ CheckDefines() { '-DHASGZ' ) export CPPTRAJ_ZLIB=$DEFINE ;; '-DHASBZ2' ) export CPPTRAJ_BZLIB=$DEFINE ;; '-DBINTRAJ' ) export CPPTRAJ_NETCDFLIB=$DEFINE ;; + '-DC11_SUPPORT' ) export CPPTRAJ_C11_SUPPORT=$DEFINE ;; '-DLIBPME' ) export CPPTRAJ_LIBPME=$DEFINE ;; '-DMPI' ) export CPPTRAJ_MPILIB=$DEFINE ;; '-DNO_MATHLIB' ) CPPTRAJ_MATHLIB='' ;; @@ -966,6 +969,7 @@ CheckEnv() { #echo "DEBUG: $DESCRIP: Checking requirement: $1" case "$1" in 'netcdf' ) TestLibrary "NetCDF" "$CPPTRAJ_NETCDFLIB" ;; + 'c++11' ) TestLibrary "C++11" "$CPPTRAJ_C11_SUPPORT" ;; 'libpme' ) TestLibrary "libPME" "$CPPTRAJ_LIBPME" ;; 'zlib' ) TestLibrary "Zlib" "$CPPTRAJ_ZLIB" ;; 'bzlib' ) TestLibrary "Bzlib" "$CPPTRAJ_BZLIB" ;; @@ -1223,6 +1227,8 @@ if [ -z "$CPPTRAJ_TEST_SETUP" ] ; then export CPPTRAJ_TEST_ERROR='Test_Error.dat' CPPTRAJ_OUTPUT='test.out' CPPTRAJ_ERROR='' + #export CPPTRAJ_RNG='--rng marsaglia' + CPPTRAJ_RNG='' # Placeholder # Process command line options CmdLineOpts $* # Determine standalone or AmberTools diff --git a/test/Test_Random/RunTest.sh b/test/Test_Random/RunTest.sh new file mode 100755 index 0000000000..25f1694790 --- /dev/null +++ b/test/Test_Random/RunTest.sh @@ -0,0 +1,54 @@ +#!/bin/bash + +. ../MasterTest.sh + +TESTNAME='Random number generator tests' + +CleanFiles rng.in random.dat mt.dat stdlib.dat pcg32.dat + +INPUT='-i rng.in' + +UNITNAME='Test Marsaglia, PCG32, and Xoshiro128++ RNGs' +cat > rng.in < rng.in < rng.in < rng.in <