Government 2013
Volume 28 Number 10A
Test Run - Implementing the National Institute of Standards and Technology Tests of Randomness Using C#
.jpg)
There’s a long history of support and cooperation between U.S. government organizations and the software development community. This article presents C# code that implements tests of randomness devised by the National Institute of Standards and Technology (NIST). In addition to providing a practical and useful addition to your programming toolset, I hope this article will persuade you that using .NET technologies to implement algorithms created by federal organizations can be significantly easier than using older technologies.
The best way to see where this article is headed is to take a look at the screenshot of a demo program in Figure 1. The demo program begins by creating and displaying a sequence of 52 binary tokens represented by 0s and 1s:
1100 0011 1010 1110 0000 1111 0000 1111 0000 1111 0000 1111 0000
.png "National Institute of Standards and Technology Tests for Randomness Demo")
Figure 1 National Institute of Standards and Technology Tests for Randomness Demo
Imagine that each token represents the color of a playing card (0 for black spades and clubs, and 1 for red hearts and diamonds), and the sequence shown is the order in which cards were dealt to you by gambling-simulation software. Does the sequence of 0s and 1s appear random, as it should be, or does the sequence look suspicious? Notice the first 16 cards seem to be rather unexceptional, but the remaining 36 cards appear to have a pattern. Did this happen by chance, or is the sequence statistically suspicious?
The C# demo program runs three tests that were devised by NIST. The first test checks the sequence for the overall frequencies of 0s and 1s. In this case, there are 25 1s and 27 0s, and the frequency test concludes that this is not suspicious. The second test examines the sequence by dividing it into six blocks of eight tokens, with the last four tokens not used. The block test also concludes there’s nothing suspicious about the pattern. The third NIST test checks the sequence for runs, which are sub-sequences in which all the tokens are the same. In this case there are 16 runs. The third test concludes that the number of runs in the input pattern is highly unlikely to have occurred if the cards were dealt randomly.
This article assumes you have at least intermediate-level skills with .NET technologies, but doesn’t assume you’re very knowledgeable about statistics. The demo program is coded in C# but you should have no trouble refactoring the demo to other .NET languages such as Visual Basic .NET or Windows PowerShell. The demo code has most error checking removed to keep the main ideas clear. The entire demo source code is available at msdn.com/magazine/msdnmaggov13.
The reference document for this article is “A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications,” NIST Special Publication 800-22, Revision 1a. I’ll refer to the document as the Test Suite. The document is available from csrc.nist.gov.
The Test Suite describes a total of 15 tests for randomness. This article presents C# implementations of the first three tests. The Test Suite was originally created before the existence of .NET technologies, so NIST provides example C language code. I’m a big fan of the C language, but C wasn’t really designed for tasks such as testing for randomness. The C download from NIST includes 10 .h files and 22 .c files, and each of the 32 files can have several functions. Even though the sample code is clearly explained, in my opinion working with a single C# file is much easier than dealing with 32 C files.
Overall Program Structure
The overall program structure of the demo, with a few minor edits, is presented in Figure 2.
Figure 2 Demo Program Structure
using System;
using System.Collections;
namespace NistRandomness
{
class NistProgram
{
static void Main(string[] args)
{
try
{
Console.WriteLine("Begin NIST tests of randomness using C# demo\n");
string bitString = "1100 0011 1010 1110 0000 1111" +
" 0000 1111 0000 1111 0000 1111 0000";
BitArray bitArray = MakeBitArray(bitString);
Console.WriteLine("Input sequence to test for randomness: \n");
ShowBitArray(bitArray, 4, 52);
Console.WriteLine("1. Testing input frequencies");
double pFreq = FrequencyTest(bitArray);
Console.WriteLine("pValue for Frequency test = " + pFreq.ToString("F4"));
if (pFreq < 0.01)
Console.WriteLine("There is evidence that sequence is NOT random");
else
Console.WriteLine("Sequence passes NIST frequency test for randomness");
int blockLength = 8;
Console.WriteLine("2. Testing input blocks (block length = " +
blockLength + ")");
double pBlock = BlockTest(bitArray, blockLength);
Console.WriteLine("pValue for Block test = " + pBlock.ToString("F4"));
if (pBlock < 0.01)
Console.WriteLine("There is evidence that sequence is NOT random");
else
Console.WriteLine("Sequence passes NIST block test for randomness");
Console.WriteLine("3. Testing input runs");
double pRuns = RunsTest(bitArray);
Console.WriteLine("pValue for Runs test = " + pRuns.ToString("F4"));
if (pRuns < 0.01)
Console.WriteLine("There is evidence that sequence is NOT random");
else
Console.WriteLine("Sequence passes NIST runs test for randomness");
Console.WriteLine("End NIST randomness demo\n");
Console.ReadLine();
}
catch (Exception ex)
{
Console.WriteLine(ex.Message);
Console.ReadLine();
}
} // Main
static BitArray MakeBitArray(string bitString) { . . }
static void ShowBitArray(BitArray bitArray, int blockSize,
int lineSize) { . . }
static double FrequencyTest(BitArray bitArray)
static double ErrorFunction(double x) { . . }
static double ErrorFunctionComplementary(double x) { . . }
static double BlockTest(BitArray bitArray, int blockLength) { . . }
static double RunsTest(BitArray bitArray) { . . }
}
public class GammaFunctions { . . }
} // ns
To create the demo, I launched Visual Studio. The demo has no significant .NET dependencies and any version of Visual Studio should work. I created a new C# console application project named NistRandomness. After the template code loaded, in the Solution Explorer window I renamed file Program.cs to the more descriptive NistProgram.cs, and Visual Studio automatically renamed class Program for me. At the top of the template-generated code, I deleted all references to namespaces except for the ones to System and Collections.
In addition to the Main method, the demo has seven static methods that implement the three tests and helper routines to store a binary sequence, plus a GammaFunction class that houses math routines used by the BlockTest method.
The demo sets up a sequence of binary tokens using a string approach:
string bitString = "1100 0011 ...";
BitArray bitArray = MakeBitArray(bitString);
In most cases, you’d be testing a sequence generated by some other system, such as a cryptographic method, or perhaps even a gambling simulation. In general, you shouldn’t have too much difficulty writing code that transforms your sequence into a string of 0s and 1s.
Calling all 15 of the NIST tests follows a consistent pattern where the input sequence is passed to a test method and the test method returns a probability value between 0 and 1. For example:
double pFreq = FrequencyTest(bitArray);
The probability represents the likelihood that the input pattern could’ve happened by chance if the pattern was in fact generated randomly. So very small values—generally those less than 0.05 or 0.01—suggest the input pattern isn’t random. The Test Suite recommends using a 0.01 significance level, but I recommend 0.05 in most cases.
The BitArray Class
The .NET BitArray class is well-suited for efficiently storing a sequence of binary tokens. Method MakeBitArray is defined in Figure 3.
Figure 3 The MakeBitArray Method
static BitArray MakeBitArray(string bitString)
{
int size = 0;
for (int i = 0; i < bitString.Length; ++i)
if (bitString[i] != ' ') ++size;
BitArray result = new BitArray(size);
int k = 0; // ptr into result
for (int i = 0; i < bitString.Length; ++i) {
if (bitString[i] == ' ') continue;
if (bitString[i] == '1') result[k] = true;
else result[k] = false;
++k;
}
return result;
}
A BitArray object is a bit unusual (no pun intended) in the sense that even though behind the scenes the contents are stored as bits, an object’s interface is exposed using Boolean true/false values. For example:
bool[] vals = new bool[] { true, false, true, false };
BitArray ba = new BitArray(vals);
This approach is extremely tedious for long input sequences. Method MakeBitArray accepts a string of 0s and 1s, and allows spaces for readability.
Method ShowBitArray accepts a BitArray object to display, and parameters to control where spaces and line breaks are placed:
static void ShowBitArray(BitArray bitArray, int blockSize,
int lineSize)
{
for (int i = 0; i < bitArray.Length; ++i) {
if (i > 0 && i % blockSize == 0) Console.Write(" ");
if (i > 0 && i % lineSize == 0) Console.WriteLine("");
if (bitArray[i] == false) Console.Write("0");
else Console.Write("1");
}
Console.WriteLine("");
}
The Frequency Test
Consider the pattern 0100 0000 0000 0010. The pattern doesn’t appear to be random because there are too many 0s. The most fundamental test for randomness is a frequency test. If a pattern is generated randomly, you’d expect the number of 0s and 1s to be roughly the same. Too many 0s or too many 1s suggest non-randomness. Method FrequencyTest computes a sum where 0s are encoded as -1 values and 1s are encoded as +1 values. If the sum is 0, there are equal numbers of 0s and 1s, but a sum different from 0, either very negative or very positive, means too many 0s or 1s. In the preceding example pattern, the sum would be -12.
The NIST Test Suite document is, in my opinion, a fantastic example of exactly how software documentation should be written: clear, concise and with excellent examples. The document explains exactly how the Frequency test works on page 2-2, and I had no problem implementing it:
static double FrequencyTest(BitArray bitArray)
{
double sum = 0;
for (int i = 0; i < bitArray.Length; ++i) {
if (bitArray[i] == false) sum = sum - 1;
else sum = sum + 1;
}
double testStat = Math.Abs(sum) / Math.Sqrt(bitArray.Length);
double rootTwo = 1.414213562373095;
double pValue = ErrorFunctionComplement(testStat / rootTwo);
return pValue;
}
The test statistic is the absolute value of the sum divided by the square root of the number of binary tokens. You can consider this a magic equation. The trickiest part of the Frequency test is the computation of the return p-value using method ErrorFunctionComplement. In math, there’s an Error Function that has a specific definition. The complement of the Error Function is just 1.0 minus the values of the Error Function. The demo defines method ErrorFunctionComplement as:
static double ErrorFunctionComplement(double x)
{
return 1 - ErrorFunction(x);
}
where:
static double ErrorFunction(double x)
{
double p = 0.3275911;
double a1 = 0.254829592;
double a2 = -0.284496736;
double a3 = 1.421413741;
double a4 = -1.453152027;
double a5 = 1.061405429;
double t = 1.0 / (1.0 + p * x);
double err = 1.0 - (((((a5 * t + a4) * t) +
a3) * t + a2) * t + a1) * t * Math.Exp(-x * x);
return err;
}
The Error Function is closely related to the Standard Normal Distribution (the bell-shaped curve) and can be used to determine how likely or unlikely a statistical event is. Method ErrorFunction uses one of many approximation equations to the true function. Here I use equation 7.1.26, which is given in the landmark reference book, “Handbook of Mathematical Functions” (Dover Publications, 1965), by Abramowitz and Stegun. The book is often just called “A&S.” This book is also a product of NIST and is indispensable if you work with numerical programming. “A&S” is freely available online at dlmf.nist.gov.
The Block Test
Consider the pattern 00000000 11111111. This pattern would pass the Frequency test because there are equal numbers of 0s and 1s, but clearly the pattern looks suspicious. The Block test is designed to address this type of non-randomness. The Block test divides a pattern into blocks and examines the number of 1s in each block. A random pattern would be expected to have about 50 percent 1s in every block. The Test Suite explains the algorithm for the Block test in detail on p. 2-4.
Briefly, the Block test accepts a block length parameter that’s the number of bits per block. From this, the number of blocks can be computed. Next, the Block test computes the proportion of 1s in each block and then uses a magic equation to compute a chi-squared test statistic. The code for method BlockTest is presented in Figure 4. An image that shows how the Block test works is shown in Figure 5.
Figure 4 The BlockTest Method
static double BlockTest(BitArray bitArray, int blockLength)
{
int numBlocks = bitArray.Length / blockLength; // 'N'
double[] proportions = new double[numBlocks];
int k = 0; // ptr into bitArray
for (int block = 0; block < numBlocks; ++block)
{
int countOnes = 0;
for (int i = 0; i < blockLength; ++i)
{
if (bitArray[k++] == true)
++countOnes;
}
proportions[block] = (countOnes * 1.0) / blockLength;
}
double summ = 0.0;
for (int block = 0; block < numBlocks; ++block)
summ = summ + (proportions[block] - 0.5) *
(proportions[block] - 0.5);
double chiSquared = 4 * blockLength * summ; // magic
double a = numBlocks / 2.0;
double x = chiSquared / 2.0;
double pValue = GammaFunctions.GammaUpper(a, x);
return pValue;
}
.png "The Block Test Algorithm")
Figure 5 The Block Test Algorithm
By far the most difficult part of the Block test—and of the entire Test Suite—is the computation of the return p-value using function GammaUpper.
Gamma Functions
It’s quite likely you’ve never encountered Gamma functions before, but they’re a cornerstone of numerical programming. There’s one primary Gamma function and several closely related functions. The Block test uses a Gamma function that’s, unfortunately, called by several different names, including the “incomplete upper Gamma function,” “GammaQ,” “upper Gamma complement” and “igamc.”
The primary Gamma function essentially extends the notion of factorial to real values. Gamma functions can be used to determine the likelihood of a complex statistical event. The demo program places five Gamma functions into a class named GammaFunctions. These functions could’ve been placed in the demo program as standalone static functions just like ErrorFunction and ErrorFunctionComplement, but housing the Gamma functions together makes them a bit more manageable. Here’s the structure of the GammaFunctions class:
public class GammaFunctions
{
public static double GammaLower(double a, double x) { . . }
public static double GammaUpper(double a, double x) { . . }
private static double LogGamma(double x) { . . }
private static double GammaLowerSer(double a, double x) { . . }
private static double GammaUpperCont(double a, double x) { . . }
}
Bear with me for a moment. The one gamma function that’s directly called by the Block test is GammaUpper. GammaUpper calls private helper functions GammaLowerSer and GammaUpperCont. Both helper functions call another helper function LogGamma. Class GammaFunctions also implements a function GammaLower that isn’t used by the Block test.
The relationships among these Gamma functions can drive you batty if you’re new to them, but you don’t need to understand these relationships to effectively use the demo code. That said, the LogGamma function computes the log of the primary Gamma function. The actual primary Gamma function is rarely computed because it overflows for even moderate input values. Function LogGamma is presented in Figure 6. The LogGamma function has been studied for decades. The approximation used in the demo employs what’s called the Lanczos algorithm.
Figure 6 The LogGamma Function
public static double LogGamma(double x)
{
double[] coef = new double[6] { 76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179E-2, -0.5395239384953E-5 };
double LogSqrtTwoPi = 0.91893853320467274178;
double denom = x + 1;
double y = x + 5.5;
double series = 1.000000000190015;
for (int i = 0; i < 6; ++i)
{
series += coef[i] / denom;
denom += 1.0;
}
return (LogSqrtTwoPi + (x + 0.5) * Math.Log(y) -
y + Math.Log(series / x));
}
Helper functions GammaLowerSer and GammaUpperCont are presented in Figure 7. The versions used here are derived from the well-known book, “Numerical Recipes in C” (Cambridge University Press, 2007), sometimes abbreviated “NR,” which in turn uses the algorithms in “A&S.” To be honest, I’m not a fan of the coding style used in “NR,” but I’ve used most of the same, rather cryptic “NR” variable names so you can exploit “NR” as a resource. The second edition of “NR” is freely available online at bit.ly/p5gRye.
Figure 7 Incomplete Gamma Function Helpers
private static double GammaLowerSer(double a, double x)
{
// Incomplete GammaLower (computed by series expansion)
if ( x < 0.0)
throw new Exception("x param less than 0.0 in GammaLowerSer");
double gln = LogGamma(a);
double ap = a;
double del = 1.0 / a;
double sum = del;
for (int n = 1; n <= 100; ++n)
{
++ap;
del *= x / ap;
sum += del;
if (Math.Abs(del) < Math.Abs(sum) * 3.0E-7) // Close enough?
return sum * Math.Exp(-x + a * Math.Log(x) - gln);
}
throw new Exception("Unable to compute GammaLowerSer " +
"to desired accuracy");
}
private static double GammaUpperCont(double a, double x)
{
// Incomplete GammaUpper computed by continuing fraction
if (x < 0.0)
throw new Exception("x param less than 0.0 in GammaUpperCont");
double gln = LogGamma(a);
double b = x + 1.0 - a;
double c = 1.0 / 1.0E-30; // Div by close to double.MinValue
double d = 1.0 / b;
double h = d;
for (int i = 1; i <= 100; ++i)
{
double an = -i * (i - a);
b += 2.0;
d = an * d + b;
if (Math.Abs(d) < 1.0E-30) d = 1.0E-30; // Got too small?
c = b + an / c;
if (Math.Abs(c) < 1.0E-30) c = 1.0E-30;
d = 1.0 / d;
double del = d * c;
h *= del;
if (Math.Abs(del - 1.0) < 3.0E-7)
return Math.Exp(-x + a * Math.Log(x) - gln) * h; // Close enough?
}
throw new Exception("Unable to compute GammaUpperCont " +
"to desired accuracy");
}
In some situations you might need to modify the code for GammaLowerSer and GammaUpperCont. Both functions are iterative approximations and stop when a maximum number of iterations (100) is reached or when an error threshold (3.0E-7 = 0.0000003) is achieved. If your Block test throws an exception, you can increase the max iterations, or increase the value of the error threshold, or both.
With the three helper functions in place, the GammaUpper function is very short:
public static double GammaUpper(double a, double x)
{
// Incomplete Gamma 'Q' (upper)
if (x < 0.0 || a <= 0.0)
throw new Exception("Bad args in GammaUpper");
if (x < a + 1)
return 1.0 - GammaLowerSer(a, x); // Indirect is faster
else
return GammaUpperCont(a, x);
}
It’s not immediately obvious what’s going on here. It turns out that GammaLowerSer and GammaUpperCont are really the same functions in the sense that the value of either of them is just 1.0 minus the value of the other. In other words, if you know one you can easily compute the other. So why two versions? As it turns out, there are two techniques to estimate an incomplete Gamma function: by a series expansion or by a continuing fraction. The series expansion is more accurate and quicker when the value of input parameter x is less than the value of input parameter a+1. The continued fraction version is more accurate and quicker otherwise.
Although it isn’t used by the Block test, the demo program implements the incomplete lower Gamma function (aka GammaP) as GammaLower:
public static double GammaLower(double a, double x)
{
// Incomplete Gamma 'P' (lower) aka 'igam'
if (x < 0.0 || a <= 0.0)
throw new Exception("Bad args in GammaLower");
if (x < a + 1)
return GammaLowerSer(a, x); // No surprise
else
return 1.0 - GammaUpperCont(a, x); // Indirectly is faster
}
The Runs Test
Consider the pattern 0101 0101 0101 0101. This pattern would pass the Frequency test because there are equal numbers of 0s and 1s. The pattern would also likely pass the Block test because each block would have roughly 50 percent of 0 bits and 50 percent of 1 bits (depending on whether the block size is even or odd). But, pretty clearly, the pattern doesn’t appear random. The Runs test catches patterns like this. A run is a sequence where consecutive bit tokens are the same. For example, the pattern 00100011 has four runs: 00, 1, 000 and 11. If a pattern is randomly generated, it’s possible to compute the expected number of runs.
The Runs test is presented in Figure 8. The test first computes the proportion of 1s in the input pattern. Then the number of runs is computed using the idea that if you walk-through a pattern, a bit at a time, every time you see a change in bit value you’ve encountered the start of a new run.
Figure 8 The Runs Test
static double RunsTest(BitArray bitArray)
{
double numOnes = 0.0;
for (int i = 0; i < bitArray.Length; ++i)
if (bitArray[i] == true)
++numOnes;
double prop = (numOnes * 1.0) / bitArray.Length;
// Double tau = 2.0 / Math.Sqrt(bitArray.Length * 1.0);
// If (Math.Abs(prop - 0.5) >= tau)
// Return 0.0; // Not-random short-circuit
int runs = 1;
for (int i = 0; i < bitArray.Length - 1; ++i)
if (bitArray[i] != bitArray[i + 1])
++runs;
double num = Math.Abs(
runs - (2 * bitArray.Length * prop * (1 - prop)));
double denom = 2 * Math.Sqrt(2.0 * bitArray.Length) *
prop * (1 - prop);
double pValue = ErrorFunctionComplement(num / denom);
return pValue;
}
The magic test statistic is computed using the overall proportion of 1s and the number of runs. The return p-value is computed using the complement of the Error Function, just as in the Frequency test.
Create Your Own .NET Tool
With the sample code and explanation in this article, you should be able to create a standalone .NET tool to analyze patterns for randomness or integrate pattern-analysis code directly into a .NET software system, using the three tests in the demo program. Additionally, you should be able to implement any or all of the other 12 tests described in the NIST Test Suite document.
Combining Strengths
In technical articles, I generally try to stick to objective information and avoid spouting subjective opinions. But I’ll step up on my soapbox this time. I’ve worked in both government public-sector software development environments (primarily projects for the U.S. Navy while living in Hawaii), and private-sector environments (primarily for Microsoft).
It’s easy for private-sector employees to take pot shots at the public sector by calling out things like massive bureaucracy, ponderous documentation and glacial development progress. And it’s easy for public-sector employees to take pot shots at the private sector by citing such factors as shoddy quality and wild-west lifecycle processes.
The tests of randomness developed by the National Institute of Standards and Technology provide a good example of how combining the relative strengths of the private and public sectors can lead to better software. I wish that every documentation resource I use was as thorough and clear as NIST documentation. And developing fairly complex mathematical software using C# is dramatically more efficient, less error-prone and quicker than using older technologies.
Dr. James McCaffrey works for Microsoft Research in Redmond, Wash. He has worked on several Microsoft products including Internet Explorer and Bing. Reach Dr. McCaffrey at jammc@microsoft.com.
Thanks to the following technical expert for reviewing this article: Dan Liebling (Microsoft Research)