From 3b9b6d0b8e7f798023c9d109c490449d528fde80 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sun, 7 Apr 2024 17:59:48 +0200 Subject: Adding upstream version 1:9.18.19. Signed-off-by: Daniel Baumann --- tests/isc/random_test.c | 787 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 787 insertions(+) create mode 100644 tests/isc/random_test.c (limited to 'tests/isc/random_test.c') 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 +#include +#include /* IWYU pragma: keep */ +#include +#include +#include +#include +#include + +#define UNIT_TESTING +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#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 -- cgit v1.2.3