random123

Counter-based Random Number Generators
git clone git://git.meso-star.com/random123.git
Log | Files | Refs | README | LICENSE

ut_Engine.cpp (8134B)


      1 /*
      2 Copyright 2010-2011, D. E. Shaw Research.
      3 All rights reserved.
      4 
      5 Redistribution and use in source and binary forms, with or without
      6 modification, are permitted provided that the following conditions are
      7 met:
      8 
      9 * Redistributions of source code must retain the above copyright
     10   notice, this list of conditions, and the following disclaimer.
     11 
     12 * Redistributions in binary form must reproduce the above copyright
     13   notice, this list of conditions, and the following disclaimer in the
     14   documentation and/or other materials provided with the distribution.
     15 
     16 * Neither the name of D. E. Shaw Research nor the names of its
     17   contributors may be used to endorse or promote products derived from
     18   this software without specific prior written permission.
     19 
     20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
     21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
     22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
     23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
     24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
     25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
     26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
     28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
     29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
     30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     31 */
     32 // TODO - really do a thorough and complete set of tests.
     33 
     34 #ifdef _MSC_FULL_VER
     35 // Engines have multiple copy constructors, quite legal C++, disable MSVC complaint
     36 #pragma warning (disable : 4521)
     37 #endif
     38 
     39 #include <Random123/philox.h>
     40 #include <Random123/aes.h>
     41 #include <Random123/threefry.h>
     42 #include <Random123/ars.h>
     43 #include <Random123/conventional/Engine.hpp>
     44 #include <Random123/ReinterpretCtr.hpp>
     45 #if R123_USE_CXX11_RANDOM
     46 #include <random>
     47 #endif
     48 #include <cassert>
     49 #include <iostream>
     50 #include <sstream>
     51 #include "util_demangle.hpp"
     52 
     53 using namespace std;
     54 using namespace r123;
     55 
     56 template <typename EType>
     57 typename EType::result_type kat1000(){
     58     // A zero return says that no KAT is known.  This makes
     59     // sense for the ReniterpretCtr-based engines which are
     60     // expected to produce endian-specific results, so we
     61     // don't have known answers for them.
     62     return 0;
     63 }
     64 
     65 #if R123_USE_64BIT
     66 #if R123_USE_PHILOX_64BIT
     67 template <> uint64_t kat1000<Engine<Philox2x64 > >(){ return R123_64BIT(10575809911605703474); }
     68 #endif
     69 template <> uint64_t kat1000<Engine<Threefry2x64 > >(){ return R123_64BIT(17578122881062615727); }
     70 #endif
     71 template <> uint32_t kat1000<Engine<Philox4x32 > >(){ return 1721865298; }
     72 template <> uint32_t kat1000<Engine<Threefry4x32 > >(){ return 874101813; }
     73 #if R123_USE_AES_OPENSSL
     74 template <> uint8_t  kat1000<Engine<AESOpenSSL16x8> >(){ return 0237; }
     75 #endif
     76 
     77 #define ASSERTEQ(A, B) assert(A==B); assert(A()==B()); assert(A()==B()); assert(A()==B())
     78 
     79 struct DummySeedSeq{
     80     template <typename ITER>
     81     void generate(ITER b, ITER e){
     82         std::fill(b, e, 1);
     83     }
     84 };
     85 
     86 template <typename EType>
     87 void doit(){
     88     EType e;
     89     cout << "doit<" << demangle(e) << ">";
     90     typedef typename EType::cbrng_type BType;
     91     typedef typename EType::result_type rtype;
     92     typedef typename BType::ctr_type ctype;
     93     typedef typename BType::key_type ktype;
     94 
     95     DummySeedSeq dummyss;
     96     EType ess(dummyss);
     97     assert(ess != e);
     98 
     99     rtype r1 = e();
    100     rtype r2 = e(); assert( r1 != r2 );
    101     rtype r3 = e(); assert( r3 != r2 && r3 != r1 );
    102 
    103     // We've elsewhere confirmed that the underlying bijections actually "work",
    104     // e.g., that they pass the Known Answer Test for some set of test vectors.
    105     // Here, we simply check that the output of the Engine corresponds, in the expected
    106     // way to the output of the underlying bijection.
    107     
    108     // Check that the first few values out of the engine correspond
    109     // to output from the underlying bijection.
    110     BType b;
    111     ctype c1 = {{}};
    112     ktype k = e.getkey();
    113     e.seed(); 
    114     for(int i=0; i<100; ++i){
    115         c1[0]++;
    116         ctype rb = b(c1, k);
    117         for(typename ctype::reverse_iterator p=rb.rbegin(); p!=rb.rend(); ++p){
    118             rtype re = e();
    119             assert( *p == re );
    120         }
    121     }
    122 
    123     // Check that discard work as expected, i.e., we can keep two
    124     // engines "in sync" by discarding from one and stepping the other.
    125     EType e2;
    126     assert(e2 != e);
    127     e2.discard(100*c1.size());
    128     ASSERTEQ(e2, e);
    129     
    130     for(int disc=1; disc<50; ++disc){
    131         e.discard(disc);
    132         assert( e != e2 );
    133         for(int j=0; j<disc; ++j) e2();
    134         ASSERTEQ(e2, e);
    135     }
    136 
    137     // Check that saving and restoring state and the copy constructor 
    138     // works as expected.
    139     ostringstream oss;
    140     oss << e2;
    141     string s2 = oss.str();
    142     int fiftyfive = 55;
    143 #if R123_USE_CXX11_TYPE_TRAITS
    144     // With CXX11, the library has type_traits to prevent
    145     // undesirable type resolution against the templated SeedSeq constructor
    146     EType e3(fiftyfive);
    147 #else
    148     // Without CXX11, we have to be careful to pass in
    149     // a bona fide rtype, and not just something that will promote
    150     // to an rtype, if we want the rtype constructor.
    151     EType e3((rtype(fiftyfive)));
    152 #endif
    153     EType esave(e);
    154     assert(e3 != e2);
    155     {
    156         istringstream iss(s2);
    157         iss >> e3;
    158     }
    159     ASSERTEQ(e3, e2);
    160     assert(e3 != esave );
    161     {
    162         istringstream iss(s2);
    163         iss >> e3;
    164     }
    165     ASSERTEQ(e3, esave);
    166     
    167     // Check that the constructor-from-rvalue works.
    168     EType e4((rtype)99);
    169     EType e5;
    170     assert(e4 != e5);
    171     assert(e4 != e3);
    172     e5.seed((rtype)99);
    173     ASSERTEQ(e4, e5);
    174 
    175 #if R123_USE_STD_RANDOM
    176     // Check that we can use an EType with a std::distribution.
    177     // Obviously, this requires <random>
    178     uniform_int_distribution<int> dieroller(1, 6);
    179     vector<int> hist(7);
    180     int NROLL = 10000;
    181     for(int i=0; i<NROLL; ++i){
    182         int roll = dieroller(e5);
    183         hist[roll]++;
    184     }
    185     double chisq = 0.;
    186     double expected = NROLL/6.;
    187     double var = NROLL*5./36.;
    188     for(int pips=1; pips<=6; ++pips){
    189         double delta = hist[pips] - expected;
    190         chisq += delta*delta/var;
    191     }
    192     // The critical value of chisq with 6 degrees of freedom
    193     // for P=0.01 is 16.81.  For P=0.05, it is 12.59
    194     const double chicrit = 12.59;
    195     if( chisq > chicrit ){
    196         printf("std::uniform_int_distribution doesn't look random.  Chisq = %g.  Does this look like the result of a fair set of dice rolls to you?\n", chisq);
    197         for(int pips=1; pips<=6; ++pips){
    198             printf("%d pips  %d times\n", pips, hist[pips]);
    199         }
    200         abort();  // a bit harsh, no?  It might just be a rare event at the 5% level...
    201     }
    202 #endif
    203 
    204     // Finally do a kat test.  
    205     EType ekat;
    206     ekat.discard(1000);
    207     typename EType::result_type r = ekat();
    208     typename EType::result_type knownanswer = kat1000<EType>();
    209     if( knownanswer != 0 && r != knownanswer )
    210         cerr << "KAT mismatch.  The 1000th random from " << demangle(ekat) << " is " << r << " it should be " << knownanswer << "\n";
    211     assert( knownanswer==0 || r == knownanswer );
    212     cout << " OK" << endl;
    213 }
    214 
    215 int main(int, char **){
    216 #if R123_USE_PHILOX_64BIT
    217     doit<Engine<Philox2x64 > >();
    218     doit<Engine<ReinterpretCtr<r123array4x32, Philox2x64 > > >();
    219 #endif
    220     doit<Engine<Philox4x32 > >();
    221     doit<Engine<Threefry4x32 > >();
    222 #if R123_USE_64BIT
    223     doit<Engine<ReinterpretCtr<r123array4x32, Threefry2x64 > > >();
    224     doit<Engine<Threefry2x64 > >();
    225     doit<Engine<ReinterpretCtr<r123array2x64, Threefry4x32 > > >();
    226 #endif
    227 
    228 #if R123_USE_AES_NI
    229     if( haveAESNI() ){
    230         doit<Engine<ARS4x32> >();
    231 #if R123_USE_64BIT
    232         doit<Engine<ReinterpretCtr<r123array2x64, ARS4x32> > >();
    233 #endif
    234         doit<Engine<AESNI4x32> >();
    235     }else{
    236         cout << "AES is compiled into the binary, but is not available on this hardware\n";
    237     }
    238 #endif
    239 #if R123_USE_AES_OPENSSL
    240     doit<Engine<AESOpenSSL16x8> >();
    241 #endif
    242 
    243     cout << "ut_Engine:  all OK" << endl;
    244     return 0;
    245 }
    246