summaryrefslogtreecommitdiffstats
path: root/tests/isc/random_test.c
diff options
context:
space:
mode:
Diffstat (limited to 'tests/isc/random_test.c')
-rw-r--r--tests/isc/random_test.c787
1 files changed, 787 insertions, 0 deletions
diff --git a/tests/isc/random_test.c b/tests/isc/random_test.c
new file mode 100644
index 0000000..1935846
--- /dev/null
+++ b/tests/isc/random_test.c
@@ -0,0 +1,787 @@
+/*
+ * Copyright (C) Internet Systems Consortium, Inc. ("ISC")
+ *
+ * SPDX-License-Identifier: MPL-2.0
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, you can obtain one at https://mozilla.org/MPL/2.0/.
+ *
+ * See the COPYRIGHT file distributed with this work for additional
+ * information regarding copyright ownership.
+ */
+
+/*
+ * IMPORTANT NOTE:
+ * These tests work by generating a large number of pseudo-random numbers
+ * and then statistically analyzing them to determine whether they seem
+ * random. The test is expected to fail on occasion by random happenstance.
+ */
+
+#include <inttypes.h>
+#include <math.h>
+#include <sched.h> /* IWYU pragma: keep */
+#include <setjmp.h>
+#include <stdarg.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+#define UNIT_TESTING
+#include <cmocka.h>
+
+#include <isc/commandline.h>
+#include <isc/mem.h>
+#include <isc/nonce.h>
+#include <isc/print.h>
+#include <isc/random.h>
+#include <isc/result.h>
+#include <isc/util.h>
+
+#include <tests/isc.h>
+
+#define REPS 25000
+
+typedef double(pvalue_func_t)(uint16_t *values, size_t length);
+
+/* igamc(), igam(), etc. were adapted (and cleaned up) from the Cephes
+ * math library:
+ *
+ * Cephes Math Library Release 2.8: June, 2000
+ * Copyright 1985, 1987, 2000 by Stephen L. Moshier
+ *
+ * The Cephes math library was released into the public domain as part
+ * of netlib.
+ */
+
+static double MACHEP = 1.11022302462515654042E-16;
+static double MAXLOG = 7.09782712893383996843E2;
+static double big = 4.503599627370496e15;
+static double biginv = 2.22044604925031308085e-16;
+
+static double
+igamc(double a, double x);
+static double
+igam(double a, double x);
+
+typedef enum {
+ ISC_RANDOM8,
+ ISC_RANDOM16,
+ ISC_RANDOM32,
+ ISC_RANDOM_BYTES,
+ ISC_RANDOM_UNIFORM,
+ ISC_NONCE_BYTES
+} isc_random_func;
+
+static double
+igamc(double a, double x) {
+ double ans, ax, c, r, t, y, z;
+ double pkm1, pkm2, qkm1, qkm2;
+
+ if ((x <= 0) || (a <= 0)) {
+ return (1.0);
+ }
+
+ if ((x < 1.0) || (x < a)) {
+ return (1.0 - igam(a, x));
+ }
+
+ ax = a * log(x) - x - lgamma(a);
+ if (ax < -MAXLOG) {
+ print_error("# igamc: UNDERFLOW, ax=%f\n", ax);
+ return (0.0);
+ }
+ ax = exp(ax);
+
+ /* continued fraction */
+ y = 1.0 - a;
+ z = x + y + 1.0;
+ c = 0.0;
+ pkm2 = 1.0;
+ qkm2 = x;
+ pkm1 = x + 1.0;
+ qkm1 = z * x;
+ ans = pkm1 / qkm1;
+
+ do {
+ double yc, pk, qk;
+ c += 1.0;
+ y += 1.0;
+ z += 2.0;
+ yc = y * c;
+ pk = pkm1 * z - pkm2 * yc;
+ qk = qkm1 * z - qkm2 * yc;
+ if (qk != 0) {
+ r = pk / qk;
+ t = fabs((ans - r) / r);
+ ans = r;
+ } else {
+ t = 1.0;
+ }
+
+ pkm2 = pkm1;
+ pkm1 = pk;
+ qkm2 = qkm1;
+ qkm1 = qk;
+
+ if (fabs(pk) > big) {
+ pkm2 *= biginv;
+ pkm1 *= biginv;
+ qkm2 *= biginv;
+ qkm1 *= biginv;
+ }
+ } while (t > MACHEP);
+
+ return (ans * ax);
+}
+
+static double
+igam(double a, double x) {
+ double ans, ax, c, r;
+
+ if ((x <= 0) || (a <= 0)) {
+ return (0.0);
+ }
+
+ if ((x > 1.0) && (x > a)) {
+ return (1.0 - igamc(a, x));
+ }
+
+ /* Compute x**a * exp(-x) / md_gamma(a) */
+ ax = a * log(x) - x - lgamma(a);
+ if (ax < -MAXLOG) {
+ print_error("# igam: UNDERFLOW, ax=%f\n", ax);
+ return (0.0);
+ }
+ ax = exp(ax);
+
+ /* power series */
+ r = a;
+ c = 1.0;
+ ans = 1.0;
+
+ do {
+ r += 1.0;
+ c *= x / r;
+ ans += c;
+ } while (c / ans > MACHEP);
+
+ return (ans * ax / a);
+}
+
+static int8_t scounts_table[65536];
+static uint8_t bitcounts_table[65536];
+
+static int8_t
+scount_calculate(uint16_t n) {
+ int i;
+ int8_t sc;
+
+ sc = 0;
+ for (i = 0; i < 16; i++) {
+ uint16_t lsb;
+
+ lsb = n & 1;
+ if (lsb != 0) {
+ sc += 1;
+ } else {
+ sc -= 1;
+ }
+
+ n >>= 1;
+ }
+
+ return (sc);
+}
+
+static uint8_t
+bitcount_calculate(uint16_t n) {
+ int i;
+ uint8_t bc;
+
+ bc = 0;
+ for (i = 0; i < 16; i++) {
+ uint16_t lsb;
+
+ lsb = n & 1;
+ if (lsb != 0) {
+ bc += 1;
+ }
+
+ n >>= 1;
+ }
+
+ return (bc);
+}
+
+static void
+tables_init(void) {
+ uint32_t i;
+
+ for (i = 0; i < 65536; i++) {
+ scounts_table[i] = scount_calculate(i);
+ bitcounts_table[i] = bitcount_calculate(i);
+ }
+}
+
+/*
+ * The following code for computing Marsaglia's rank is based on the
+ * implementation in cdbinrnk.c from the diehard tests by George
+ * Marsaglia.
+ *
+ * This function destroys (modifies) the data passed in bits.
+ */
+static uint32_t
+matrix_binaryrank(uint32_t *bits, size_t rows, size_t cols) {
+ unsigned int rt = 0;
+ uint32_t rank = 0;
+ uint32_t tmp;
+
+ for (size_t k = 0; k < rows; k++) {
+ size_t i = k;
+
+ while (rt >= cols || ((bits[i] >> rt) & 1) == 0) {
+ i++;
+
+ if (i < rows) {
+ continue;
+ } else {
+ rt++;
+ if (rt < cols) {
+ i = k;
+ continue;
+ }
+ }
+
+ return (rank);
+ }
+
+ rank++;
+ if (i != k) {
+ tmp = bits[i];
+ bits[i] = bits[k];
+ bits[k] = tmp;
+ }
+
+ for (size_t j = i + 1; j < rows; j++) {
+ if (((bits[j] >> rt) & 1) == 0) {
+ continue;
+ } else {
+ bits[j] ^= bits[k];
+ }
+ }
+
+ rt++;
+ }
+
+ return (rank);
+}
+
+static void
+random_test(pvalue_func_t *func, isc_random_func test_func) {
+ uint32_t m;
+ uint32_t j;
+ uint32_t histogram[11] = { 0 };
+ uint32_t passed;
+ double proportion;
+ double p_hat;
+ double lower_confidence, higher_confidence;
+ double chi_square;
+ double p_value_t;
+ double alpha;
+
+ tables_init();
+
+ m = 1000;
+ passed = 0;
+
+ for (j = 0; j < m; j++) {
+ uint32_t i;
+ uint32_t values[REPS];
+ uint16_t *uniform_values;
+ double p_value;
+
+ switch (test_func) {
+ case ISC_RANDOM8:
+ for (i = 0; i < (sizeof(values) / sizeof(*values)); i++)
+ {
+ values[i] = isc_random8();
+ }
+ break;
+ case ISC_RANDOM16:
+ for (i = 0; i < (sizeof(values) / sizeof(*values)); i++)
+ {
+ values[i] = isc_random16();
+ }
+ break;
+ case ISC_RANDOM32:
+ for (i = 0; i < (sizeof(values) / sizeof(*values)); i++)
+ {
+ values[i] = isc_random32();
+ }
+ break;
+ case ISC_RANDOM_BYTES:
+ isc_random_buf(values, sizeof(values));
+ break;
+ case ISC_RANDOM_UNIFORM:
+ uniform_values = (uint16_t *)values;
+ for (i = 0;
+ i < (sizeof(values) / (sizeof(*uniform_values)));
+ i++)
+ {
+ uniform_values[i] =
+ isc_random_uniform(UINT16_MAX);
+ }
+ break;
+ case ISC_NONCE_BYTES:
+ isc_nonce_buf(values, sizeof(values));
+ break;
+ }
+
+ p_value = (*func)((uint16_t *)values, REPS * 2);
+ if (p_value >= 0.01) {
+ passed++;
+ }
+
+ assert_in_range(p_value, 0.0, 1.0);
+
+ i = (int)floor(p_value * 10);
+ histogram[i]++;
+ }
+
+ /*
+ * Check proportion of sequences passing a test (see section
+ * 4.2.1 in NIST SP 800-22).
+ */
+ alpha = 0.01; /* the significance level */
+ proportion = (double)passed / (double)m;
+ p_hat = 1.0 - alpha;
+ lower_confidence = p_hat - (3.0 * sqrt((p_hat * (1.0 - p_hat)) / m));
+ higher_confidence = p_hat + (3.0 * sqrt((p_hat * (1.0 - p_hat)) / m));
+
+ assert_in_range(proportion, lower_confidence, higher_confidence);
+
+ /*
+ * Check uniform distribution of p-values (see section 4.2.2 in
+ * NIST SP 800-22).
+ */
+
+ /* Fold histogram[10] (p_value = 1.0) into histogram[9] for
+ * interval [0.9, 1.0]
+ */
+ histogram[9] += histogram[10];
+ histogram[10] = 0;
+
+ /* Pre-requisite that at least 55 sequences are processed. */
+ assert_true(m >= 55);
+
+ chi_square = 0.0;
+ for (j = 0; j < 10; j++) {
+ double numer;
+ double denom;
+
+ numer = (histogram[j] - (m / 10.0)) *
+ (histogram[j] - (m / 10.0));
+ denom = m / 10.0;
+ chi_square += numer / denom;
+ }
+
+ p_value_t = igamc(9 / 2.0, chi_square / 2.0);
+
+ assert_true(p_value_t >= 0.0001);
+}
+
+/*
+ * This is a frequency (monobits) test taken from the NIST SP 800-22
+ * RANDOM test suite.
+ */
+static double
+monobit(uint16_t *values, size_t length) {
+ size_t i;
+ int32_t scount;
+ uint32_t numbits;
+ double s_obs;
+ double p_value;
+
+ UNUSED(mctx);
+
+ numbits = length * sizeof(*values) * 8;
+ scount = 0;
+
+ for (i = 0; i < length; i++) {
+ scount += scounts_table[values[i]];
+ }
+
+ /* Preconditions (section 2.1.7 in NIST SP 800-22) */
+ assert_true(numbits >= 100);
+
+ s_obs = abs(scount) / sqrt(numbits);
+ p_value = erfc(s_obs / sqrt(2.0));
+
+ return (p_value);
+}
+
+/*
+ * This is the runs test taken from the NIST SP 800-22 RNG test suite.
+ */
+static double
+runs(uint16_t *values, size_t length) {
+ size_t i;
+ uint32_t bcount;
+ uint32_t numbits;
+ double pi;
+ double tau;
+ uint32_t j;
+ uint32_t b;
+ uint8_t bit_prev;
+ uint32_t v_obs;
+ double numer;
+ double denom;
+ double p_value;
+
+ UNUSED(mctx);
+
+ numbits = length * sizeof(*values) * 8;
+ bcount = 0;
+
+ for (i = 0; i < length; i++) {
+ bcount += bitcounts_table[values[i]];
+ }
+
+ pi = (double)bcount / (double)numbits;
+ tau = 2.0 / sqrt(numbits);
+
+ /* Preconditions (section 2.3.7 in NIST SP 800-22) */
+ assert_true(numbits >= 100);
+
+ /*
+ * Pre-condition implied from the monobit test. This can fail
+ * for some sequences, and the p-value is taken as 0 in these
+ * cases.
+ */
+ if (fabs(pi - 0.5) >= tau) {
+ return (0.0);
+ }
+
+ /* Compute v_obs */
+ j = 0;
+ b = 14;
+ bit_prev = (values[j] & (1U << 15)) == 0 ? 0 : 1;
+
+ v_obs = 0;
+
+ for (i = 1; i < numbits; i++) {
+ uint8_t bit_this = (values[j] & (1U << b)) == 0 ? 0 : 1;
+ if (b == 0) {
+ b = 15;
+ j++;
+ } else {
+ b--;
+ }
+
+ v_obs += bit_this ^ bit_prev;
+
+ bit_prev = bit_this;
+ }
+
+ v_obs += 1;
+
+ numer = fabs(v_obs - (2.0 * numbits * pi * (1.0 - pi)));
+ denom = 2.0 * sqrt(2.0 * numbits) * pi * (1.0 - pi);
+
+ p_value = erfc(numer / denom);
+
+ return (p_value);
+}
+
+/*
+ * This is the block frequency test taken from the NIST SP 800-22 RNG
+ * test suite.
+ */
+static double
+blockfrequency(uint16_t *values, size_t length) {
+ uint32_t i;
+ uint32_t numbits;
+ uint32_t mbits;
+ uint32_t mwords;
+ uint32_t numblocks;
+ double *pi;
+ double chi_square;
+ double p_value;
+
+ numbits = length * sizeof(*values) * 8;
+ mbits = 32000;
+ mwords = mbits / 16;
+ numblocks = numbits / mbits;
+
+ /* Preconditions (section 2.2.7 in NIST SP 800-22) */
+ assert_true(numbits >= 100);
+ assert_true(mbits >= 20);
+ assert_true((double)mbits > (0.01 * numbits));
+ assert_true(numblocks < 100);
+ assert_true(numbits >= (mbits * numblocks));
+
+ pi = isc_mem_get(mctx, numblocks * sizeof(double));
+ assert_non_null(pi);
+
+ for (i = 0; i < numblocks; i++) {
+ uint32_t j;
+ pi[i] = 0.0;
+ for (j = 0; j < mwords; j++) {
+ uint32_t idx;
+
+ idx = i * mwords + j;
+ pi[i] += bitcounts_table[values[idx]];
+ }
+ pi[i] /= mbits;
+ }
+
+ /* Compute chi_square */
+ chi_square = 0.0;
+ for (i = 0; i < numblocks; i++) {
+ chi_square += (pi[i] - 0.5) * (pi[i] - 0.5);
+ }
+
+ chi_square *= 4 * mbits;
+
+ isc_mem_put(mctx, pi, numblocks * sizeof(double));
+
+ p_value = igamc(numblocks * 0.5, chi_square * 0.5);
+
+ return (p_value);
+}
+
+/*
+ * This is the binary matrix rank test taken from the NIST SP 800-22 RNG
+ * test suite.
+ */
+static double
+binarymatrixrank(uint16_t *values, size_t length) {
+ uint32_t i;
+ size_t matrix_m;
+ size_t matrix_q;
+ uint32_t num_matrices;
+ size_t numbits;
+ uint32_t fm_0;
+ uint32_t fm_1;
+ uint32_t fm_rest;
+ double term1;
+ double term2;
+ double term3;
+ double chi_square;
+ double p_value;
+
+ UNUSED(mctx);
+
+ matrix_m = 32;
+ matrix_q = 32;
+ num_matrices = length / ((matrix_m * matrix_q) / 16);
+ numbits = num_matrices * matrix_m * matrix_q;
+
+ /* Preconditions (section 2.5.7 in NIST SP 800-22) */
+ assert_int_equal(matrix_m, 32);
+ assert_int_equal(matrix_q, 32);
+ assert_true(numbits >= (38 * matrix_m * matrix_q));
+
+ fm_0 = 0;
+ fm_1 = 0;
+ fm_rest = 0;
+ for (i = 0; i < num_matrices; i++) {
+ /*
+ * Each uint32_t supplies 32 bits, so a 32x32 bit matrix
+ * takes up uint32_t array of size 32.
+ */
+ uint32_t bits[32];
+ int j;
+ uint32_t rank;
+
+ for (j = 0; j < 32; j++) {
+ size_t idx;
+ uint32_t r1;
+ uint32_t r2;
+
+ idx = i * ((matrix_m * matrix_q) / 16);
+ idx += j * 2;
+
+ r1 = values[idx];
+ r2 = values[idx + 1];
+ bits[j] = (r1 << 16) | r2;
+ }
+
+ rank = matrix_binaryrank(bits, matrix_m, matrix_q);
+
+ if (rank == matrix_m) {
+ fm_0++;
+ } else if (rank == (matrix_m - 1)) {
+ fm_1++;
+ } else {
+ fm_rest++;
+ }
+ }
+
+ /* Compute chi_square */
+ term1 = ((fm_0 - (0.2888 * num_matrices)) *
+ (fm_0 - (0.2888 * num_matrices))) /
+ (0.2888 * num_matrices);
+ term2 = ((fm_1 - (0.5776 * num_matrices)) *
+ (fm_1 - (0.5776 * num_matrices))) /
+ (0.5776 * num_matrices);
+ term3 = ((fm_rest - (0.1336 * num_matrices)) *
+ (fm_rest - (0.1336 * num_matrices))) /
+ (0.1336 * num_matrices);
+
+ chi_square = term1 + term2 + term3;
+
+ p_value = exp(-chi_square * 0.5);
+
+ return (p_value);
+}
+
+/***
+ *** Tests for isc_random32() function
+ ***/
+
+/* Monobit test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random32_monobit) {
+ UNUSED(state);
+
+ random_test(monobit, ISC_RANDOM32);
+}
+
+/* Runs test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random32_runs) {
+ UNUSED(state);
+
+ random_test(runs, ISC_RANDOM32);
+}
+
+/* Block frequency test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random32_blockfrequency) {
+ UNUSED(state);
+
+ random_test(blockfrequency, ISC_RANDOM32);
+}
+
+/* Binary matrix rank test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random32_binarymatrixrank) {
+ UNUSED(state);
+
+ random_test(binarymatrixrank, ISC_RANDOM32);
+}
+
+/***
+ *** Tests for isc_random_bytes() function
+ ***/
+
+/* Monobit test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_bytes_monobit) {
+ UNUSED(state);
+
+ random_test(monobit, ISC_RANDOM_BYTES);
+}
+
+/* Runs test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_bytes_runs) {
+ UNUSED(state);
+
+ random_test(runs, ISC_RANDOM_BYTES);
+}
+
+/* Block frequency test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_bytes_blockfrequency) {
+ UNUSED(state);
+
+ random_test(blockfrequency, ISC_RANDOM_BYTES);
+}
+
+/* Binary matrix rank test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_bytes_binarymatrixrank) {
+ UNUSED(state);
+
+ random_test(binarymatrixrank, ISC_RANDOM_BYTES);
+}
+
+/***
+ *** Tests for isc_random_uniform() function:
+ ***/
+
+/* Monobit test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_uniform_monobit) {
+ UNUSED(state);
+
+ random_test(monobit, ISC_RANDOM_UNIFORM);
+}
+
+/* Runs test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_uniform_runs) {
+ UNUSED(state);
+
+ random_test(runs, ISC_RANDOM_UNIFORM);
+}
+
+/* Block frequency test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_uniform_blockfrequency) {
+ UNUSED(state);
+
+ random_test(blockfrequency, ISC_RANDOM_UNIFORM);
+}
+
+/* Binary matrix rank test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_random_uniform_binarymatrixrank) {
+ UNUSED(state);
+
+ random_test(binarymatrixrank, ISC_RANDOM_UNIFORM);
+}
+
+/* Tests for isc_nonce_bytes() function */
+
+/* Monobit test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_nonce_bytes_monobit) {
+ UNUSED(state);
+
+ random_test(monobit, ISC_NONCE_BYTES);
+}
+
+/* Runs test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_nonce_bytes_runs) {
+ UNUSED(state);
+
+ random_test(runs, ISC_NONCE_BYTES);
+}
+
+/* Block frequency test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_nonce_bytes_blockfrequency) {
+ UNUSED(state);
+
+ random_test(blockfrequency, ISC_NONCE_BYTES);
+}
+
+/* Binary matrix rank test for the RANDOM */
+ISC_RUN_TEST_IMPL(isc_nonce_bytes_binarymatrixrank) {
+ UNUSED(state);
+
+ random_test(binarymatrixrank, ISC_NONCE_BYTES);
+}
+
+ISC_TEST_LIST_START
+
+ISC_TEST_ENTRY(isc_random32_monobit)
+ISC_TEST_ENTRY(isc_random32_runs)
+ISC_TEST_ENTRY(isc_random32_blockfrequency)
+ISC_TEST_ENTRY(isc_random32_binarymatrixrank)
+ISC_TEST_ENTRY(isc_random_bytes_monobit)
+ISC_TEST_ENTRY(isc_random_bytes_runs)
+ISC_TEST_ENTRY(isc_random_bytes_blockfrequency)
+ISC_TEST_ENTRY(isc_random_bytes_binarymatrixrank)
+ISC_TEST_ENTRY(isc_random_uniform_monobit)
+ISC_TEST_ENTRY(isc_random_uniform_runs)
+ISC_TEST_ENTRY(isc_random_uniform_blockfrequency)
+ISC_TEST_ENTRY(isc_random_uniform_binarymatrixrank)
+ISC_TEST_ENTRY(isc_nonce_bytes_monobit)
+ISC_TEST_ENTRY(isc_nonce_bytes_runs)
+ISC_TEST_ENTRY(isc_nonce_bytes_blockfrequency)
+ISC_TEST_ENTRY(isc_nonce_bytes_binarymatrixrank)
+
+ISC_TEST_LIST_END
+
+ISC_TEST_MAIN
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-03 20:19:25 +0000