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