How to Test with TestU01

In an earlier post I showed you how to test a PRNG with PractRand. Now it's the turn of the venerable TestU01. It's a little more finicky, but quite doable. Let's see…

Downloading and Building TestU01

You can find the TestU01 website here. It configures in the traditional way for Unix software.

I'm going to assume you don't want to do a full install of TestU01 as a general library on your system, you just want to use it briefly, so we'll do a quick and easy temporary install. The commands below are all you need to do to build a working test setup with TestU01 on a Unix system (Linux, macOS, or the Windows subsystem for Linux),

mkdir TestU01
cd TestU01
basedir=`pwd`
curl -OL http://simul.iro.umontreal.ca/testu01/TestU01.zip
unzip -q TestU01.zip
cd TestU01-1.2.3
./configure --prefix="$basedir"
make -j 6
make -j 6 install
cd ..
mkdir lib-so
mv lib/*.so lib-so/.

The last two lines aren't 100% necessary, but dealing with dynamic libraries is more trouble than it's worth here. This way we'll use the static libraries.

Writing a Test Program to test Xorshift-32

The TestU01 documentation explains how to test a simple 32-bit generator, and gives sample code.

// Adapted from TestU01 manual, Figure 2.2, Figure 2.4

#include "TestU01.h"

// Example PRNG: Xorshift 32

static unsigned int y = 2463534242U;

unsigned int xorshift (void)
{
    y ^= (y << 13);
    y ^= (y >> 17);
    return y ^= (y << 5);
}

int main()
{
    // Create TestU01 PRNG object for our generator
    unif01_Gen* gen = unif01_CreateExternGenBits("Xorshift 32", xorshift);

    // Run the tests.
    bbattery_SmallCrush(gen);

    // Clean up.
    unif01_DeleteExternGenBits(gen);

    return 0;
}

To compile it, we need to run

gcc -std=c99 -Wall -O3 -o test-xorshift-minimal test-xorshift-minimal.c -Iinclude -Llib -ltestu01 -lprobdist -lmylib -lm

It's vital that we include the -Iinclude -Llib -ltestu01 -lprobdist -lmylib -lm part to let the compiler find the TestU01 libraries and link with them.

If we run the resulting executable, we get a lot of output as it goes through each test, culminating in this summary:

========= Summary results of SmallCrush =========

 Version:          TestU01 1.2.3
 Generator:        Xorshift 32
 Number of statistics:  15
 Total CPU time:   00:00:05.46
 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  BirthdaySpacings                 eps
  2  Collision                      1 - eps1
  6  MaxOft                         6.7e-16
  8  MatrixRank                       eps
 10  RandomWalk1 H                   5.7e-7
 ----------------------------------------------
 All other tests were passed

Less than six seconds of testing, and we see that unimproved XorShift generators have problems. (Truncated XorShift*, on the other hand, would have passed with flying colors.)

As short and sweet as this test code was, it's not ideal in a few ways. We can't run the code with a different seed, we can't easily change what test we run, and it doesn't also test the reversed bits (which is important with TestU01). We'll address these issues in the next examples.

Testing the 32-bit Mersenne Twister (C++)

When you add some command line options to allow the user to specify a seed, and the ability to choose test options, and so on, the code balloons a bit, but it's nothing earth shattering in its complexity:

#include <random>
#include <string>
#include <cstdio>
#include <cstdint>
#include <cstddef>

extern "C" {
    #include "TestU01.h"
}

// GENERATOR SECTION
//
// This is the only part of the code that needs to change to switch to
// a new generator.

const char* gen_name = "std mt19937";  // TestU01 doesn't like colons!!?!

const int MAX_SEEDS = 1;
uint64_t seed_data[MAX_SEEDS];

uint32_t gen32()
{
    static std::mt19937 rng(seed_data[0]);

    return rng();
}

// END OF GENERATOR SECTION

inline uint32_t rev32(uint32_t v)
{
    // https://graphics.stanford.edu/~seander/bithacks.html
    // swap odd and even bits
    v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
    // swap consecutive pairs
    v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
    // swap nibbles ...
    v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
    // swap bytes
    v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
    // swap 2-byte-long pairs
    v = ( v >> 16             ) | ( v               << 16);
    return v;
}

uint32_t gen32_rev()
{
    return rev32(gen32());
}

const char* progname;

void usage()
{
    printf("%s: [-v] [-r] [seeds...]\n", progname);
    exit(1);
}

int main (int argc, char** argv)
{
    progname = argv[0];

    // Config options for TestU01
    swrite_Basic = FALSE;  // Turn of TestU01 verbosity by default
                           // reenable by -v option.

    // Config options for generator output
    bool reverseBits = false;

    // Config options for tests
    bool testSmallCrush = false;
    bool testCrush = false;
    bool testBigCrush = false;
    bool testLinComp = false;

    // Handle command-line option switches
    while (1) {
        --argc; ++argv;
        if ((argc == 0) || (argv[0][0] != '-'))
            break;
        if ((argv[0][1]=='\0') || (argv[0][2]!='\0'))
            usage();
        switch(argv[0][1]) {
        case 'r':
            reverseBits = true;
            break;      
        case 's':
            testSmallCrush = true;
            break;
        case 'm':
            testCrush = true;
            break;
        case 'b':
            testBigCrush = true;
            break;
        case 'l':
            testLinComp = true;
            break;
        case 'v':
            swrite_Basic = TRUE;
            break;
        default:
            usage();
        }
    }

    // Name of the generator

    std::string genName = gen_name;
    if (reverseBits)
        genName += " [Reversed]";

    // Determine a default test if need be

    if (!(testSmallCrush || testCrush || testBigCrush || testLinComp)) {
        testCrush = true;
    }

    // Initialize seed-data array, either using command-line arguments
    // or std::random_device.

    printf("Testing %s:\n", genName.c_str());
    printf("- seed_data[%u] = { ", MAX_SEEDS);
    std::random_device rdev;
    for (int i = 0; i < MAX_SEEDS; ++i) {
        if (argc >= i+1) {
            seed_data[i] = strtoull(argv[i],0,0);
        } else {
            seed_data[i] = (uint64_t(rdev()) << 32) | rdev();
        }
        printf("%s0x%016lx", i == 0 ? "" : ", ", seed_data[i]);
    }
    printf("}\n");
    fflush(stdout);

    // Create a generator for TestU01.

    unif01_Gen* gen =
        unif01_CreateExternGenBits((char*) genName.c_str(),
                                   reverseBits ? gen32 : gen32_rev);

    // Run tests.

    if (testSmallCrush) {
        bbattery_SmallCrush(gen);
        fflush(stdout);
    }
    if (testCrush) {
        bbattery_Crush(gen);
        fflush(stdout);
    }
    if (testBigCrush) {
        bbattery_BigCrush(gen);
        fflush(stdout);
    }
    if (testLinComp) {
        scomp_Res* res = scomp_CreateRes();
        swrite_Basic = TRUE;
        for (int size : {250, 500, 1000, 5000, 25000, 50000, 75000})
            scomp_LinearComp(gen, res, 1, size, 0, 1);
        scomp_DeleteRes(res);
        fflush(stdout);
    }

    // Clean up.

    unif01_DeleteExternGenBits(gen);

    return 0;
}

We can compile this code with

g++ -std=c++14 -Wall -O3 -o test-mt19937 test-mt19937.cpp -Iinclude -Llib -ltestu01 -lprobdist -lmylib -lm

And run it with, for example,

linux$ ./test-mt19937 -r -s
Testing std mt19937 [Reversed]:
- seed_data[1] = { 0xc37415b2ce04157e}

========= Summary results of SmallCrush =========

 Version:          TestU01 1.2.3
 Generator:        std mt19937 [Reversed]
 Number of statistics:  15
 Total CPU time:   00:00:06.37

 All tests were passed

One neat thing about this code is that by running it with the -l option, we can reveal the same “Linear Complexity” test failure that the BigCrush and Crush tests reveal, but do so in a few seconds rather than having to wait hours for the same result. We'll show output from this test when we test the 64-bit Mersenne Twister.

Testing the 64-bit Mersenne Twister (C++)

Testing 64-bit generators with TestU01 is a bit of a pain, because we can only test 32 bits at a time. But we can adapt the above code into a version for 64-bit generators that handles the issue.

#include <random>
#include <string>
#include <cstdio>
#include <cstdint>
#include <cstddef>

extern "C" {
    #include "TestU01.h"
}

// GENERATOR SECTION
//
// This is the only part of the code that needs to change to switch to
// a new generator.

const char* gen_name = "std mt19937";  // TestU01 doesn't like colons!!?!

const int MAX_SEEDS = 1;
uint64_t seed_data[MAX_SEEDS];

uint32_t gen32()
{
    static std::mt19937 rng(seed_data[0]);

    return rng();
}

// END OF GENERATOR SECTION

inline uint32_t rev32(uint32_t v)
{
    // https://graphics.stanford.edu/~seander/bithacks.html
    // swap odd and even bits
    v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
    // swap consecutive pairs
    v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
    // swap nibbles ...
    v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
    // swap bytes
    v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
    // swap 2-byte-long pairs
    v = ( v >> 16             ) | ( v               << 16);
    return v;
}

uint32_t gen32_rev()
{
    return rev32(gen32());
}

const char* progname;

void usage()
{
    printf("%s: [-v] [-r] [seeds...]\n", progname);
    exit(1);
}

int main (int argc, char** argv)
{
    progname = argv[0];

    // Config options for TestU01
    swrite_Basic = FALSE;  // Turn of TestU01 verbosity by default
                           // reenable by -v option.

    // Config options for generator output
    bool reverseBits = false;

    // Config options for tests
    bool testSmallCrush = false;
    bool testCrush = false;
    bool testBigCrush = false;
    bool testLinComp = false;

    // Handle command-line option switches
    while (1) {
        --argc; ++argv;
        if ((argc == 0) || (argv[0][0] != '-'))
            break;
        if ((argv[0][1]=='\0') || (argv[0][2]!='\0'))
            usage();
        switch(argv[0][1]) {
        case 'r':
            reverseBits = true;
            break;      
        case 's':
            testSmallCrush = true;
            break;
        case 'm':
            testCrush = true;
            break;
        case 'b':
            testBigCrush = true;
            break;
        case 'l':
            testLinComp = true;
            break;
        case 'v':
            swrite_Basic = TRUE;
            break;
        default:
            usage();
        }
    }

    // Name of the generator

    std::string genName = gen_name;
    if (reverseBits)
        genName += " [Reversed]";

    // Determine a default test if need be

    if (!(testSmallCrush || testCrush || testBigCrush || testLinComp)) {
        testCrush = true;
    }

    // Initialize seed-data array, either using command-line arguments
    // or std::random_device.

    printf("Testing %s:\n", genName.c_str());
    printf("- seed_data[%u] = { ", MAX_SEEDS);
    std::random_device rdev;
    for (int i = 0; i < MAX_SEEDS; ++i) {
        if (argc >= i+1) {
            seed_data[i] = strtoull(argv[i],0,0);
        } else {
            seed_data[i] = (uint64_t(rdev()) << 32) | rdev();
        }
        printf("%s0x%016lx", i == 0 ? "" : ", ", seed_data[i]);
    }
    printf("}\n");
    fflush(stdout);

    // Create a generator for TestU01.

    unif01_Gen* gen =
        unif01_CreateExternGenBits((char*) genName.c_str(),
                                   reverseBits ? gen32 : gen32_rev);

    // Run tests.

    if (testSmallCrush) {
        bbattery_SmallCrush(gen);
        fflush(stdout);
    }
    if (testCrush) {
        bbattery_Crush(gen);
        fflush(stdout);
    }
    if (testBigCrush) {
        bbattery_BigCrush(gen);
        fflush(stdout);
    }
    if (testLinComp) {
        scomp_Res* res = scomp_CreateRes();
        swrite_Basic = TRUE;
        for (int size : {250, 500, 1000, 5000, 25000, 50000, 75000})
            scomp_LinearComp(gen, res, 1, size, 0, 1);
        scomp_DeleteRes(res);
        fflush(stdout);
    }

    // Clean up.

    unif01_DeleteExternGenBits(gen);

    return 0;
}

Again, we compile with

g++ -std=c++14 -Wall -O3 -o test-mt19937_64 test-mt19937_64.cpp -Iinclude -Llib -ltestu01 -lprobdist -lmylib -lm

and here's a run (piped through a couple of greps to strip blank and unimportant lines):

$ ./test-mt19937_64 -h -l | grep . | egrep -v '^(HOST|Generator)'
Testing std mt19937_64 [High bits]:
- seed_data[1] = { 0x1981b0ad2c8efea9}
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 250,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    1
Chi2 statistic for size of jumps      :    1.17
p-value of test                       :    0.28
-----------------------------------------------
Normal statistic for number of jumps  :    1.11
p-value of test                       :    0.13
-----------------------------------------------
CPU time used                    :  00:00:00.00
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 500,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    2
Chi2 statistic for size of jumps      :    2.20
p-value of test                       :    0.33
-----------------------------------------------
Normal statistic for number of jumps  :    1.22
p-value of test                       :    0.11
-----------------------------------------------
CPU time used                    :  00:00:00.00
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 1000,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    3
Chi2 statistic for size of jumps      :    3.57
p-value of test                       :    0.31
-----------------------------------------------
Normal statistic for number of jumps  :    0.95
p-value of test                       :    0.17
-----------------------------------------------
CPU time used                    :  00:00:00.00
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 5000,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    5
Chi2 statistic for size of jumps      :    1.78
p-value of test                       :    0.88
-----------------------------------------------
Normal statistic for number of jumps  :  -0.053
p-value of test                       :    0.52
-----------------------------------------------
CPU time used                    :  00:00:00.01
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 25000,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    8
Chi2 statistic for size of jumps      :    9.45
p-value of test                       :    0.31
-----------------------------------------------
Normal statistic for number of jumps  :    1.44
p-value of test                       :    0.07
-----------------------------------------------
CPU time used                    :  00:00:00.45
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 50000,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    9
Chi2 statistic for size of jumps      :    5.63
p-value of test                       :    0.78
-----------------------------------------------
Normal statistic for number of jumps  :  -32.26
p-value of test                       : 1 - eps1    *****
-----------------------------------------------
CPU time used                    :  00:00:01.35
***********************************************************
std mt19937_64 [High bits]
scomp_LinearComp test:
-----------------------------------------------
   N =  1,  n = 75000,  r =  0,    s = 1
-----------------------------------------------
Number of degrees of freedom          :    9
Chi2 statistic for size of jumps      :   10.25
p-value of test                       :    0.33
-----------------------------------------------
Normal statistic for number of jumps  :  -91.51
p-value of test                       : 1 - eps1    *****
-----------------------------------------------
CPU time used                    :  00:00:01.68

Here, in 1.68 seconds of testing, we learn exactly how/where the Mersenne Twister fails. It survives up to 50,000 outputs, but then fails. We could run BigCrush instead, but we'd wait more than three hours to learn less about this failure. If there were other failures, we'd learn about them, too, but there aren't any; this is the only one. Linear complexity is the test that catches out all (G)LFSR generators that have insufficiently strong output functions.

Conclusion

Using a random-number test suite isn't hard. You can do it, too!