1 /* $NetBSD: random_test.c,v 1.3 2025/01/26 16:25:50 christos Exp $ */ 2 3 /* 4 * Copyright (C) Internet Systems Consortium, Inc. ("ISC") 5 * 6 * SPDX-License-Identifier: MPL-2.0 7 * 8 * This Source Code Form is subject to the terms of the Mozilla Public 9 * License, v. 2.0. If a copy of the MPL was not distributed with this 10 * file, you can obtain one at https://mozilla.org/MPL/2.0/. 11 * 12 * See the COPYRIGHT file distributed with this work for additional 13 * information regarding copyright ownership. 14 */ 15 16 /* 17 * IMPORTANT NOTE: 18 * These tests work by generating a large number of pseudo-random numbers 19 * and then statistically analyzing them to determine whether they seem 20 * random. The test is expected to fail on occasion by random happenstance. 21 */ 22 23 #include <inttypes.h> 24 #include <math.h> 25 #include <sched.h> /* IWYU pragma: keep */ 26 #include <setjmp.h> 27 #include <stdarg.h> 28 #include <stddef.h> 29 #include <stdlib.h> 30 #include <string.h> 31 32 #define UNIT_TESTING 33 #include <cmocka.h> 34 35 #include <isc/commandline.h> 36 #include <isc/mem.h> 37 #include <isc/nonce.h> 38 #include <isc/random.h> 39 #include <isc/result.h> 40 #include <isc/util.h> 41 42 #include <tests/isc.h> 43 44 #define REPS 25000 45 46 typedef double(pvalue_func_t)(uint16_t *values, size_t length); 47 48 /* igamc(), igam(), etc. were adapted (and cleaned up) from the Cephes 49 * math library: 50 * 51 * Cephes Math Library Release 2.8: June, 2000 52 * Copyright 1985, 1987, 2000 by Stephen L. Moshier 53 * 54 * The Cephes math library was released into the public domain as part 55 * of netlib. 56 */ 57 58 static double MACHEP = 1.11022302462515654042E-16; 59 static double MAXLOG = 7.09782712893383996843E2; 60 static double big = 4.503599627370496e15; 61 static double biginv = 2.22044604925031308085e-16; 62 63 static double 64 igamc(double a, double x); 65 static double 66 igam(double a, double x); 67 68 typedef enum { 69 ISC_RANDOM8, 70 ISC_RANDOM16, 71 ISC_RANDOM32, 72 ISC_RANDOM_BYTES, 73 ISC_RANDOM_UNIFORM, 74 ISC_NONCE_BYTES 75 } isc_random_func; 76 77 static double 78 igamc(double a, double x) { 79 double ans, ax, c, r, t, y, z; 80 double pkm1, pkm2, qkm1, qkm2; 81 82 if ((x <= 0) || (a <= 0)) { 83 return 1.0; 84 } 85 86 if ((x < 1.0) || (x < a)) { 87 return 1.0 - igam(a, x); 88 } 89 90 ax = a * log(x) - x - lgamma(a); 91 if (ax < -MAXLOG) { 92 print_error("# igamc: UNDERFLOW, ax=%f\n", ax); 93 return 0.0; 94 } 95 ax = exp(ax); 96 97 /* continued fraction */ 98 y = 1.0 - a; 99 z = x + y + 1.0; 100 c = 0.0; 101 pkm2 = 1.0; 102 qkm2 = x; 103 pkm1 = x + 1.0; 104 qkm1 = z * x; 105 ans = pkm1 / qkm1; 106 107 do { 108 double yc, pk, qk; 109 c += 1.0; 110 y += 1.0; 111 z += 2.0; 112 yc = y * c; 113 pk = pkm1 * z - pkm2 * yc; 114 qk = qkm1 * z - qkm2 * yc; 115 if (qk != 0) { 116 r = pk / qk; 117 t = fabs((ans - r) / r); 118 ans = r; 119 } else { 120 t = 1.0; 121 } 122 123 pkm2 = pkm1; 124 pkm1 = pk; 125 qkm2 = qkm1; 126 qkm1 = qk; 127 128 if (fabs(pk) > big) { 129 pkm2 *= biginv; 130 pkm1 *= biginv; 131 qkm2 *= biginv; 132 qkm1 *= biginv; 133 } 134 } while (t > MACHEP); 135 136 return ans * ax; 137 } 138 139 static double 140 igam(double a, double x) { 141 double ans, ax, c, r; 142 143 if ((x <= 0) || (a <= 0)) { 144 return 0.0; 145 } 146 147 if ((x > 1.0) && (x > a)) { 148 return 1.0 - igamc(a, x); 149 } 150 151 /* Compute x**a * exp(-x) / md_gamma(a) */ 152 ax = a * log(x) - x - lgamma(a); 153 if (ax < -MAXLOG) { 154 print_error("# igam: UNDERFLOW, ax=%f\n", ax); 155 return 0.0; 156 } 157 ax = exp(ax); 158 159 /* power series */ 160 r = a; 161 c = 1.0; 162 ans = 1.0; 163 164 do { 165 r += 1.0; 166 c *= x / r; 167 ans += c; 168 } while (c / ans > MACHEP); 169 170 return ans * ax / a; 171 } 172 173 static int8_t scounts_table[65536]; 174 static uint8_t bitcounts_table[65536]; 175 176 static int8_t 177 scount_calculate(uint16_t n) { 178 int i; 179 int8_t sc; 180 181 sc = 0; 182 for (i = 0; i < 16; i++) { 183 uint16_t lsb; 184 185 lsb = n & 1; 186 if (lsb != 0) { 187 sc += 1; 188 } else { 189 sc -= 1; 190 } 191 192 n >>= 1; 193 } 194 195 return sc; 196 } 197 198 static uint8_t 199 bitcount_calculate(uint16_t n) { 200 int i; 201 uint8_t bc; 202 203 bc = 0; 204 for (i = 0; i < 16; i++) { 205 uint16_t lsb; 206 207 lsb = n & 1; 208 if (lsb != 0) { 209 bc += 1; 210 } 211 212 n >>= 1; 213 } 214 215 return bc; 216 } 217 218 static void 219 tables_init(void) { 220 uint32_t i; 221 222 for (i = 0; i < 65536; i++) { 223 scounts_table[i] = scount_calculate(i); 224 bitcounts_table[i] = bitcount_calculate(i); 225 } 226 } 227 228 /* 229 * The following code for computing Marsaglia's rank is based on the 230 * implementation in cdbinrnk.c from the diehard tests by George 231 * Marsaglia. 232 * 233 * This function destroys (modifies) the data passed in bits. 234 */ 235 static uint32_t 236 matrix_binaryrank(uint32_t *bits, size_t rows, size_t cols) { 237 unsigned int rt = 0; 238 uint32_t rank = 0; 239 uint32_t tmp; 240 241 for (size_t k = 0; k < rows; k++) { 242 size_t i = k; 243 244 while (rt >= cols || ((bits[i] >> rt) & 1) == 0) { 245 i++; 246 247 if (i < rows) { 248 continue; 249 } else { 250 rt++; 251 if (rt < cols) { 252 i = k; 253 continue; 254 } 255 } 256 257 return rank; 258 } 259 260 rank++; 261 if (i != k) { 262 tmp = bits[i]; 263 bits[i] = bits[k]; 264 bits[k] = tmp; 265 } 266 267 for (size_t j = i + 1; j < rows; j++) { 268 if (((bits[j] >> rt) & 1) == 0) { 269 continue; 270 } else { 271 bits[j] ^= bits[k]; 272 } 273 } 274 275 rt++; 276 } 277 278 return rank; 279 } 280 281 static void 282 random_test(pvalue_func_t *func, isc_random_func test_func) { 283 uint32_t m; 284 uint32_t j; 285 uint32_t histogram[11] = { 0 }; 286 uint32_t passed; 287 double proportion; 288 double p_hat; 289 double lower_confidence, higher_confidence; 290 double chi_square; 291 double p_value_t; 292 double alpha; 293 294 tables_init(); 295 296 m = 1000; 297 passed = 0; 298 299 for (j = 0; j < m; j++) { 300 uint32_t i; 301 uint32_t values[REPS]; 302 uint16_t *uniform_values; 303 double p_value; 304 305 switch (test_func) { 306 case ISC_RANDOM8: 307 for (i = 0; i < (sizeof(values) / sizeof(*values)); i++) 308 { 309 values[i] = isc_random8(); 310 } 311 break; 312 case ISC_RANDOM16: 313 for (i = 0; i < (sizeof(values) / sizeof(*values)); i++) 314 { 315 values[i] = isc_random16(); 316 } 317 break; 318 case ISC_RANDOM32: 319 for (i = 0; i < (sizeof(values) / sizeof(*values)); i++) 320 { 321 values[i] = isc_random32(); 322 } 323 break; 324 case ISC_RANDOM_BYTES: 325 isc_random_buf(values, sizeof(values)); 326 break; 327 case ISC_RANDOM_UNIFORM: 328 uniform_values = (uint16_t *)values; 329 for (i = 0; 330 i < (sizeof(values) / (sizeof(*uniform_values))); 331 i++) 332 { 333 uniform_values[i] = 334 isc_random_uniform(UINT16_MAX); 335 } 336 break; 337 case ISC_NONCE_BYTES: 338 isc_nonce_buf(values, sizeof(values)); 339 break; 340 } 341 342 p_value = (*func)((uint16_t *)values, REPS * 2); 343 if (p_value >= 0.01) { 344 passed++; 345 } 346 347 assert_in_range(p_value, 0.0, 1.0); 348 349 i = (int)floor(p_value * 10); 350 histogram[i]++; 351 } 352 353 /* 354 * Check proportion of sequences passing a test (see section 355 * 4.2.1 in NIST SP 800-22). 356 */ 357 alpha = 0.01; /* the significance level */ 358 proportion = (double)passed / (double)m; 359 p_hat = 1.0 - alpha; 360 lower_confidence = p_hat - (3.0 * sqrt((p_hat * (1.0 - p_hat)) / m)); 361 higher_confidence = p_hat + (3.0 * sqrt((p_hat * (1.0 - p_hat)) / m)); 362 363 assert_in_range(proportion, lower_confidence, higher_confidence); 364 365 /* 366 * Check uniform distribution of p-values (see section 4.2.2 in 367 * NIST SP 800-22). 368 */ 369 370 /* Fold histogram[10] (p_value = 1.0) into histogram[9] for 371 * interval [0.9, 1.0] 372 */ 373 histogram[9] += histogram[10]; 374 histogram[10] = 0; 375 376 /* Pre-requisite that at least 55 sequences are processed. */ 377 assert_true(m >= 55); 378 379 chi_square = 0.0; 380 for (j = 0; j < 10; j++) { 381 double numer; 382 double denom; 383 384 numer = (histogram[j] - (m / 10.0)) * 385 (histogram[j] - (m / 10.0)); 386 denom = m / 10.0; 387 chi_square += numer / denom; 388 } 389 390 p_value_t = igamc(9 / 2.0, chi_square / 2.0); 391 392 assert_true(p_value_t >= 0.0001); 393 } 394 395 /* 396 * This is a frequency (monobits) test taken from the NIST SP 800-22 397 * RANDOM test suite. 398 */ 399 static double 400 monobit(uint16_t *values, size_t length) { 401 size_t i; 402 int32_t scount; 403 uint32_t numbits; 404 double s_obs; 405 double p_value; 406 407 UNUSED(mctx); 408 409 numbits = length * sizeof(*values) * 8; 410 scount = 0; 411 412 for (i = 0; i < length; i++) { 413 scount += scounts_table[values[i]]; 414 } 415 416 /* Preconditions (section 2.1.7 in NIST SP 800-22) */ 417 assert_true(numbits >= 100); 418 419 s_obs = abs(scount) / sqrt(numbits); 420 p_value = erfc(s_obs / sqrt(2.0)); 421 422 return p_value; 423 } 424 425 /* 426 * This is the runs test taken from the NIST SP 800-22 RNG test suite. 427 */ 428 static double 429 runs(uint16_t *values, size_t length) { 430 size_t i; 431 uint32_t bcount; 432 uint32_t numbits; 433 double pi; 434 double tau; 435 uint32_t j; 436 uint32_t b; 437 uint8_t bit_prev; 438 uint32_t v_obs; 439 double numer; 440 double denom; 441 double p_value; 442 443 UNUSED(mctx); 444 445 numbits = length * sizeof(*values) * 8; 446 bcount = 0; 447 448 for (i = 0; i < length; i++) { 449 bcount += bitcounts_table[values[i]]; 450 } 451 452 pi = (double)bcount / (double)numbits; 453 tau = 2.0 / sqrt(numbits); 454 455 /* Preconditions (section 2.3.7 in NIST SP 800-22) */ 456 assert_true(numbits >= 100); 457 458 /* 459 * Pre-condition implied from the monobit test. This can fail 460 * for some sequences, and the p-value is taken as 0 in these 461 * cases. 462 */ 463 if (fabs(pi - 0.5) >= tau) { 464 return 0.0; 465 } 466 467 /* Compute v_obs */ 468 j = 0; 469 b = 14; 470 bit_prev = (values[j] & (1U << 15)) == 0 ? 0 : 1; 471 472 v_obs = 0; 473 474 for (i = 1; i < numbits; i++) { 475 uint8_t bit_this = (values[j] & (1U << b)) == 0 ? 0 : 1; 476 if (b == 0) { 477 b = 15; 478 j++; 479 } else { 480 b--; 481 } 482 483 v_obs += bit_this ^ bit_prev; 484 485 bit_prev = bit_this; 486 } 487 488 v_obs += 1; 489 490 numer = fabs(v_obs - (2.0 * numbits * pi * (1.0 - pi))); 491 denom = 2.0 * sqrt(2.0 * numbits) * pi * (1.0 - pi); 492 493 p_value = erfc(numer / denom); 494 495 return p_value; 496 } 497 498 /* 499 * This is the block frequency test taken from the NIST SP 800-22 RNG 500 * test suite. 501 */ 502 static double 503 blockfrequency(uint16_t *values, size_t length) { 504 uint32_t i; 505 uint32_t numbits; 506 uint32_t mbits; 507 uint32_t mwords; 508 uint32_t numblocks; 509 double *pi; 510 double chi_square; 511 double p_value; 512 513 numbits = length * sizeof(*values) * 8; 514 mbits = 32000; 515 mwords = mbits / 16; 516 numblocks = numbits / mbits; 517 518 /* Preconditions (section 2.2.7 in NIST SP 800-22) */ 519 assert_true(numbits >= 100); 520 assert_true(mbits >= 20); 521 assert_true((double)mbits > (0.01 * numbits)); 522 assert_true(numblocks < 100); 523 assert_true(numbits >= (mbits * numblocks)); 524 525 pi = isc_mem_cget(mctx, numblocks, sizeof(double)); 526 assert_non_null(pi); 527 528 for (i = 0; i < numblocks; i++) { 529 uint32_t j; 530 pi[i] = 0.0; 531 for (j = 0; j < mwords; j++) { 532 uint32_t idx; 533 534 idx = i * mwords + j; 535 pi[i] += bitcounts_table[values[idx]]; 536 } 537 pi[i] /= mbits; 538 } 539 540 /* Compute chi_square */ 541 chi_square = 0.0; 542 for (i = 0; i < numblocks; i++) { 543 chi_square += (pi[i] - 0.5) * (pi[i] - 0.5); 544 } 545 546 chi_square *= 4 * mbits; 547 548 isc_mem_cput(mctx, pi, numblocks, sizeof(double)); 549 550 p_value = igamc(numblocks * 0.5, chi_square * 0.5); 551 552 return p_value; 553 } 554 555 /* 556 * This is the binary matrix rank test taken from the NIST SP 800-22 RNG 557 * test suite. 558 */ 559 static double 560 binarymatrixrank(uint16_t *values, size_t length) { 561 uint32_t i; 562 size_t matrix_m; 563 size_t matrix_q; 564 uint32_t num_matrices; 565 size_t numbits; 566 uint32_t fm_0; 567 uint32_t fm_1; 568 uint32_t fm_rest; 569 double term1; 570 double term2; 571 double term3; 572 double chi_square; 573 double p_value; 574 575 UNUSED(mctx); 576 577 matrix_m = 32; 578 matrix_q = 32; 579 num_matrices = length / ((matrix_m * matrix_q) / 16); 580 numbits = num_matrices * matrix_m * matrix_q; 581 582 /* Preconditions (section 2.5.7 in NIST SP 800-22) */ 583 assert_int_equal(matrix_m, 32); 584 assert_int_equal(matrix_q, 32); 585 assert_true(numbits >= (38 * matrix_m * matrix_q)); 586 587 fm_0 = 0; 588 fm_1 = 0; 589 fm_rest = 0; 590 for (i = 0; i < num_matrices; i++) { 591 /* 592 * Each uint32_t supplies 32 bits, so a 32x32 bit matrix 593 * takes up uint32_t array of size 32. 594 */ 595 uint32_t bits[32]; 596 int j; 597 uint32_t rank; 598 599 for (j = 0; j < 32; j++) { 600 size_t idx; 601 uint32_t r1; 602 uint32_t r2; 603 604 idx = i * ((matrix_m * matrix_q) / 16); 605 idx += j * 2; 606 607 r1 = values[idx]; 608 r2 = values[idx + 1]; 609 bits[j] = (r1 << 16) | r2; 610 } 611 612 rank = matrix_binaryrank(bits, matrix_m, matrix_q); 613 614 if (rank == matrix_m) { 615 fm_0++; 616 } else if (rank == (matrix_m - 1)) { 617 fm_1++; 618 } else { 619 fm_rest++; 620 } 621 } 622 623 /* Compute chi_square */ 624 term1 = ((fm_0 - (0.2888 * num_matrices)) * 625 (fm_0 - (0.2888 * num_matrices))) / 626 (0.2888 * num_matrices); 627 term2 = ((fm_1 - (0.5776 * num_matrices)) * 628 (fm_1 - (0.5776 * num_matrices))) / 629 (0.5776 * num_matrices); 630 term3 = ((fm_rest - (0.1336 * num_matrices)) * 631 (fm_rest - (0.1336 * num_matrices))) / 632 (0.1336 * num_matrices); 633 634 chi_square = term1 + term2 + term3; 635 636 p_value = exp(-chi_square * 0.5); 637 638 return p_value; 639 } 640 641 /*** 642 *** Tests for isc_random32() function 643 ***/ 644 645 /* Ensure the RNG has been automatically seeded. */ 646 ISC_RUN_TEST_IMPL(isc_random32_initialized) { 647 UNUSED(state); 648 649 assert_int_not_equal(isc_random32(), 0); 650 } 651 652 /* Monobit test for the RANDOM */ 653 ISC_RUN_TEST_IMPL(isc_random32_monobit) { 654 UNUSED(state); 655 656 random_test(monobit, ISC_RANDOM32); 657 } 658 659 /* Runs test for the RANDOM */ 660 ISC_RUN_TEST_IMPL(isc_random32_runs) { 661 UNUSED(state); 662 663 random_test(runs, ISC_RANDOM32); 664 } 665 666 /* Block frequency test for the RANDOM */ 667 ISC_RUN_TEST_IMPL(isc_random32_blockfrequency) { 668 UNUSED(state); 669 670 random_test(blockfrequency, ISC_RANDOM32); 671 } 672 673 /* Binary matrix rank test for the RANDOM */ 674 ISC_RUN_TEST_IMPL(isc_random32_binarymatrixrank) { 675 UNUSED(state); 676 677 random_test(binarymatrixrank, ISC_RANDOM32); 678 } 679 680 /*** 681 *** Tests for isc_random_bytes() function 682 ***/ 683 684 /* Monobit test for the RANDOM */ 685 ISC_RUN_TEST_IMPL(isc_random_bytes_monobit) { 686 UNUSED(state); 687 688 random_test(monobit, ISC_RANDOM_BYTES); 689 } 690 691 /* Runs test for the RANDOM */ 692 ISC_RUN_TEST_IMPL(isc_random_bytes_runs) { 693 UNUSED(state); 694 695 random_test(runs, ISC_RANDOM_BYTES); 696 } 697 698 /* Block frequency test for the RANDOM */ 699 ISC_RUN_TEST_IMPL(isc_random_bytes_blockfrequency) { 700 UNUSED(state); 701 702 random_test(blockfrequency, ISC_RANDOM_BYTES); 703 } 704 705 /* Binary matrix rank test for the RANDOM */ 706 ISC_RUN_TEST_IMPL(isc_random_bytes_binarymatrixrank) { 707 UNUSED(state); 708 709 random_test(binarymatrixrank, ISC_RANDOM_BYTES); 710 } 711 712 /*** 713 *** Tests for isc_random_uniform() function: 714 ***/ 715 716 /* Monobit test for the RANDOM */ 717 ISC_RUN_TEST_IMPL(isc_random_uniform_monobit) { 718 UNUSED(state); 719 720 random_test(monobit, ISC_RANDOM_UNIFORM); 721 } 722 723 /* Runs test for the RANDOM */ 724 ISC_RUN_TEST_IMPL(isc_random_uniform_runs) { 725 UNUSED(state); 726 727 random_test(runs, ISC_RANDOM_UNIFORM); 728 } 729 730 /* Block frequency test for the RANDOM */ 731 ISC_RUN_TEST_IMPL(isc_random_uniform_blockfrequency) { 732 UNUSED(state); 733 734 random_test(blockfrequency, ISC_RANDOM_UNIFORM); 735 } 736 737 /* Binary matrix rank test for the RANDOM */ 738 ISC_RUN_TEST_IMPL(isc_random_uniform_binarymatrixrank) { 739 UNUSED(state); 740 741 random_test(binarymatrixrank, ISC_RANDOM_UNIFORM); 742 } 743 744 /* Tests for isc_nonce_bytes() function */ 745 746 /* Monobit test for the RANDOM */ 747 ISC_RUN_TEST_IMPL(isc_nonce_bytes_monobit) { 748 UNUSED(state); 749 750 random_test(monobit, ISC_NONCE_BYTES); 751 } 752 753 /* Runs test for the RANDOM */ 754 ISC_RUN_TEST_IMPL(isc_nonce_bytes_runs) { 755 UNUSED(state); 756 757 random_test(runs, ISC_NONCE_BYTES); 758 } 759 760 /* Block frequency test for the RANDOM */ 761 ISC_RUN_TEST_IMPL(isc_nonce_bytes_blockfrequency) { 762 UNUSED(state); 763 764 random_test(blockfrequency, ISC_NONCE_BYTES); 765 } 766 767 /* Binary matrix rank test for the RANDOM */ 768 ISC_RUN_TEST_IMPL(isc_nonce_bytes_binarymatrixrank) { 769 UNUSED(state); 770 771 random_test(binarymatrixrank, ISC_NONCE_BYTES); 772 } 773 774 ISC_TEST_LIST_START 775 776 ISC_TEST_ENTRY(isc_random32_initialized) 777 ISC_TEST_ENTRY(isc_random32_monobit) 778 ISC_TEST_ENTRY(isc_random32_runs) 779 ISC_TEST_ENTRY(isc_random32_blockfrequency) 780 ISC_TEST_ENTRY(isc_random32_binarymatrixrank) 781 ISC_TEST_ENTRY(isc_random_bytes_monobit) 782 ISC_TEST_ENTRY(isc_random_bytes_runs) 783 ISC_TEST_ENTRY(isc_random_bytes_blockfrequency) 784 ISC_TEST_ENTRY(isc_random_bytes_binarymatrixrank) 785 ISC_TEST_ENTRY(isc_random_uniform_monobit) 786 ISC_TEST_ENTRY(isc_random_uniform_runs) 787 ISC_TEST_ENTRY(isc_random_uniform_blockfrequency) 788 ISC_TEST_ENTRY(isc_random_uniform_binarymatrixrank) 789 ISC_TEST_ENTRY(isc_nonce_bytes_monobit) 790 ISC_TEST_ENTRY(isc_nonce_bytes_runs) 791 ISC_TEST_ENTRY(isc_nonce_bytes_blockfrequency) 792 ISC_TEST_ENTRY(isc_nonce_bytes_binarymatrixrank) 793 794 ISC_TEST_LIST_END 795 796 ISC_TEST_MAIN 797