random_test.c

Go to the documentation of this file.
00001 /*
00002  * Copyright (C) 2014  Internet Systems Consortium, Inc. ("ISC")
00003  *
00004  * Permission to use, copy, modify, and/or distribute this software for any
00005  * purpose with or without fee is hereby granted, provided that the above
00006  * copyright notice and this permission notice appear in all copies.
00007  *
00008  * THE SOFTWARE IS PROVIDED "AS IS" AND ISC DISCLAIMS ALL WARRANTIES WITH
00009  * REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
00010  * AND FITNESS.  IN NO EVENT SHALL ISC BE LIABLE FOR ANY SPECIAL, DIRECT,
00011  * INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
00012  * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
00013  * OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
00014  * PERFORMANCE OF THIS SOFTWARE.
00015  */
00016 
00017 #include <config.h>
00018 
00019 #include <isc/random.h>
00020 #include <isc/result.h>
00021 #include <isc/mem.h>
00022 #include <isc/util.h>
00023 
00024 #include <atf-c.h>
00025 
00026 #include <stdlib.h>
00027 #include <stdint.h>
00028 #include <math.h>
00029 
00030 typedef double (pvalue_func_t)(isc_mem_t *mctx,
00031                                isc_uint16_t *values, size_t length);
00032 
00033 /* igamc(), igam(), etc. were adapted (and cleaned up) from the Cephes
00034  * math library:
00035  *
00036  * Cephes Math Library Release 2.8:  June, 2000
00037  * Copyright 1985, 1987, 2000 by Stephen L. Moshier
00038  *
00039  * The Cephes math library was released into the public domain as part
00040  * of netlib.
00041 */
00042 
00043 static double MACHEP = 1.11022302462515654042E-16;
00044 static double MAXLOG = 7.09782712893383996843E2;
00045 static double big = 4.503599627370496e15;
00046 static double biginv =  2.22044604925031308085e-16;
00047 
00048 static double igamc(double a, double x);
00049 static double igam(double a, double x);
00050 
00051 static double
00052 igamc(double a, double x) {
00053         double ans, ax, c, yc, r, t, y, z;
00054         double pk, pkm1, pkm2, qk, qkm1, qkm2;
00055 
00056         if ((x <= 0) || (a <= 0))
00057                 return (1.0);
00058 
00059         if ((x < 1.0) || (x < a))
00060                 return (1.0 - igam(a, x));
00061 
00062         ax = a * log(x) - x - lgamma(a);
00063         if (ax < -MAXLOG) {
00064                 fprintf(stderr, "igamc: UNDERFLOW, ax=%f\n", ax);
00065                 return (0.0);
00066         }
00067         ax = exp(ax);
00068 
00069         /* continued fraction */
00070         y = 1.0 - a;
00071         z = x + y + 1.0;
00072         c = 0.0;
00073         pkm2 = 1.0;
00074         qkm2 = x;
00075         pkm1 = x + 1.0;
00076         qkm1 = z * x;
00077         ans = pkm1 / qkm1;
00078 
00079         do {
00080                 c += 1.0;
00081                 y += 1.0;
00082                 z += 2.0;
00083                 yc = y * c;
00084                 pk = pkm1 * z  -  pkm2 * yc;
00085                 qk = qkm1 * z  -  qkm2 * yc;
00086                 if (qk != 0) {
00087                         r = pk / qk;
00088                         t = fabs((ans - r) / r);
00089                         ans = r;
00090                 } else
00091                         t = 1.0;
00092 
00093                 pkm2 = pkm1;
00094                 pkm1 = pk;
00095                 qkm2 = qkm1;
00096                 qkm1 = qk;
00097 
00098                 if (fabs(pk) > big) {
00099                         pkm2 *= biginv;
00100                         pkm1 *= biginv;
00101                         qkm2 *= biginv;
00102                         qkm1 *= biginv;
00103                 }
00104         } while (t > MACHEP);
00105 
00106         return (ans * ax);
00107 }
00108 
00109 static double
00110 igam(double a, double x) {
00111         double ans, ax, c, r;
00112 
00113         if ((x <= 0) || (a <= 0))
00114                 return (0.0);
00115 
00116         if ((x > 1.0) && (x > a))
00117                 return (1.0 - igamc(a, x));
00118 
00119         /* Compute  x**a * exp(-x) / md_gamma(a)  */
00120         ax = a * log(x) - x - lgamma(a);
00121         if( ax < -MAXLOG ) {
00122                 fprintf(stderr, "igam: UNDERFLOW, ax=%f\n", ax);
00123                 return (0.0);
00124         }
00125         ax = exp(ax);
00126 
00127         /* power series */
00128         r = a;
00129         c = 1.0;
00130         ans = 1.0;
00131 
00132         do {
00133                 r += 1.0;
00134                 c *= x / r;
00135                 ans += c;
00136         } while (c / ans > MACHEP);
00137 
00138         return (ans * ax / a);
00139 }
00140 
00141 static isc_int8_t scounts_table[65536];
00142 static isc_uint8_t bitcounts_table[65536];
00143 
00144 static isc_int8_t
00145 scount_calculate(isc_uint16_t n) {
00146         int i;
00147         isc_int8_t sc;
00148 
00149         sc = 0;
00150         for (i = 0; i < 16; i++) {
00151                 isc_uint16_t lsb;
00152 
00153                 lsb = n & 1;
00154                 if (lsb != 0)
00155                         sc += 1;
00156                 else
00157                         sc -= 1;
00158 
00159                 n >>= 1;
00160         }
00161 
00162         return (sc);
00163 }
00164 
00165 static isc_uint8_t
00166 bitcount_calculate(isc_uint16_t n) {
00167         int i;
00168         isc_uint8_t bc;
00169 
00170         bc = 0;
00171         for (i = 0; i < 16; i++) {
00172                 isc_uint16_t lsb;
00173 
00174                 lsb = n & 1;
00175                 if (lsb != 0)
00176                         bc += 1;
00177 
00178                 n >>= 1;
00179         }
00180 
00181         return (bc);
00182 }
00183 
00184 static void
00185 tables_init(void) {
00186         isc_uint32_t i;
00187 
00188         for (i = 0; i < 65536; i++) {
00189                 scounts_table[i] = scount_calculate(i);
00190                 bitcounts_table[i] = bitcount_calculate(i);
00191         }
00192 }
00193 
00194 /*
00195  * The following code for computing Marsaglia's rank is based on the
00196  * implementation in cdbinrnk.c from the diehard tests by George
00197  * Marsaglia.
00198  *
00199  * This function destroys (modifies) the data passed in bits.
00200  */
00201 static isc_uint32_t
00202 matrix_binaryrank(isc_uint32_t *bits, ssize_t rows, ssize_t cols) {
00203         ssize_t i, j, k;
00204         int rt = 0;
00205         isc_uint32_t rank = 0;
00206         isc_uint32_t tmp;
00207 
00208         for (k = 0; k < rows; k++) {
00209                 i = k;
00210 
00211                 while (((bits[i] >> rt) & 1) == 0) {
00212                         i++;
00213 
00214                         if (i < rows)
00215                                 continue;
00216                         else {
00217                                 rt++;
00218                                 if (rt < cols) {
00219                                         i = k;
00220                                         continue;
00221                                 }
00222                         }
00223 
00224                         return (rank);
00225                 }
00226 
00227                 rank++;
00228                 if (i != k) {
00229                         tmp = bits[i];
00230                         bits[i] = bits[k];
00231                         bits[k] = tmp;
00232                 }
00233 
00234                 for (j = i + 1; j < rows; j++) {
00235                         if (((bits[j] >> rt) & 1) == 0)
00236                                 continue;
00237                         else
00238                                 bits[j] ^= bits[k];
00239                 }
00240 
00241                 rt++;
00242         }
00243 
00244         return (rank);
00245 }
00246 
00247 static void
00248 random_test(pvalue_func_t *func) {
00249         isc_mem_t *mctx = NULL;
00250         isc_result_t result;
00251         isc_rng_t *rng;
00252         isc_uint32_t m;
00253         isc_uint32_t j;
00254         isc_uint32_t histogram[11];
00255         isc_uint32_t passed;
00256         double proportion;
00257         double p_hat;
00258         double lower_confidence;
00259         double chi_square;
00260         double p_value_t;
00261 
00262         tables_init();
00263 
00264         result = isc_mem_create(0, 0, &mctx);
00265         ATF_REQUIRE_EQ(result, ISC_R_SUCCESS);
00266 
00267         rng = NULL;
00268         result = isc_rng_create(mctx, NULL, &rng);
00269         ATF_REQUIRE_EQ(result, ISC_R_SUCCESS);
00270 
00271         m = 1000;
00272         passed = 0;
00273 
00274         for (j = 0; j < 11; j++)
00275                 histogram[j] = 0;
00276 
00277         for (j = 0; j < m; j++) {
00278                 isc_uint32_t i;
00279                 isc_uint16_t values[128000];
00280                 double p_value;
00281 
00282                 for (i = 0; i < 128000; i++)
00283                         values[i] = isc_rng_random(rng);
00284 
00285                 p_value = (*func)(mctx, values, 128000);
00286                 if (p_value >= 0.01)
00287                         passed++;
00288 
00289                 ATF_REQUIRE(p_value >= 0.0);
00290                 ATF_REQUIRE(p_value <= 1.0);
00291 
00292                 i = (int) floor(p_value * 10);
00293                 histogram[i]++;
00294         }
00295 
00296         isc_rng_detach(&rng);
00297 
00298         /* Fold histogram[10] (p_value = 1.0) into histogram[9] for
00299          * interval [0.9, 1.0]
00300          */
00301         histogram[9] += histogram[10];
00302         histogram[10] = 0;
00303 
00304         /*
00305          * Check proportion of sequences passing a test (see section
00306          * 4.2.1 in NIST SP 800-22).
00307          */
00308         proportion = (double) passed / (double) m;
00309         p_hat = 1 - 0.01; /* alpha is 0.01 in the NIST tests */
00310         lower_confidence = p_hat - (3.0 * sqrt((p_hat * (1 - p_hat)) / m));
00311 
00312         /* Debug message, not displayed when running via atf-run */
00313         printf("passed=%u/1000\n", passed);
00314         printf("lower_confidence=%f, proportion=%f\n",
00315                lower_confidence, proportion);
00316 
00317         ATF_REQUIRE(proportion >= lower_confidence);
00318 
00319         /*
00320          * Check uniform distribution of p-values (see section 4.2.2 in
00321          * NIST SP 800-22).
00322          */
00323 
00324         /* Pre-requisite that at least 55 sequences are processed. */
00325         ATF_REQUIRE(m >= 55);
00326 
00327         chi_square = 0.0;
00328         for (j = 0; j < 10; j++) {
00329                 double numer;
00330                 double denom;
00331 
00332                 /* Debug message, not displayed when running via atf-run */
00333                 printf("hist%u=%u ", j, histogram[j]);
00334 
00335                 numer = (histogram[j] - (m / 10.0)) *
00336                         (histogram[j] - (m / 10.0));
00337                 denom = m / 10.0;
00338                 chi_square += numer / denom;
00339         }
00340 
00341         printf("\n");
00342 
00343         p_value_t = igamc(9 / 2.0, chi_square / 2.0);
00344 
00345         ATF_REQUIRE(p_value_t >= 0.0001);
00346 }
00347 
00348 /*
00349  * This is a frequency (monobits) test taken from the NIST SP 800-22
00350  * RNG test suite.
00351  */
00352 static double
00353 monobit(isc_mem_t *mctx, isc_uint16_t *values, size_t length) {
00354         size_t i;
00355         isc_int32_t scount;
00356         isc_uint32_t numbits;
00357         double s_obs;
00358         double p_value;
00359 
00360         UNUSED(mctx);
00361 
00362         numbits = length * 16;
00363         scount = 0;
00364 
00365         for (i = 0; i < length; i++)
00366                 scount += scounts_table[values[i]];
00367 
00368         /* Preconditions (section 2.1.7 in NIST SP 800-22) */
00369         ATF_REQUIRE(numbits >= 100);
00370 
00371         /* Debug message, not displayed when running via atf-run */
00372         printf("numbits=%u, scount=%d\n", numbits, scount);
00373 
00374         s_obs = abs(scount) / sqrt(numbits);
00375         p_value = erfc(s_obs / sqrt(2.0));
00376 
00377         return (p_value);
00378 }
00379 
00380 /*
00381  * This is the runs test taken from the NIST SP 800-22 RNG test suite.
00382  */
00383 static double
00384 runs(isc_mem_t *mctx, isc_uint16_t *values, size_t length) {
00385         size_t i;
00386         isc_uint32_t bcount;
00387         isc_uint32_t numbits;
00388         double pi;
00389         double tau;
00390         isc_uint32_t j;
00391         isc_uint32_t b;
00392         isc_uint8_t bit_this;
00393         isc_uint8_t bit_prev;
00394         isc_uint32_t v_obs;
00395         double numer;
00396         double denom;
00397         double p_value;
00398 
00399         UNUSED(mctx);
00400 
00401         numbits = length * 16;
00402         bcount = 0;
00403 
00404         for (i = 0; i < 128000; i++)
00405                 bcount += bitcounts_table[values[i]];
00406 
00407         /* Debug message, not displayed when running via atf-run */
00408         printf("numbits=%u, bcount=%u\n", numbits, bcount);
00409 
00410         pi = (double) bcount / (double) numbits;
00411         tau = 2.0 / sqrt(numbits);
00412 
00413         /* Preconditions (section 2.3.7 in NIST SP 800-22) */
00414         ATF_REQUIRE(numbits >= 100);
00415 
00416         /*
00417          * Pre-condition implied from the monobit test. This can fail
00418          * for some sequences, and the p-value is taken as 0 in these
00419          * cases.
00420          */
00421         if (fabs(pi - 0.5) >= tau)
00422                 return (0.0);
00423 
00424         /* Compute v_obs */
00425         j = 0;
00426         b = 14;
00427         bit_prev = (values[j] & (1U << 15)) == 0 ? 0 : 1;
00428 
00429         v_obs = 0;
00430 
00431         for (i = 1; i < numbits; i++) {
00432                 bit_this = (values[j] & (1U << b)) == 0 ? 0 : 1;
00433                 if (b == 0) {
00434                         b = 15;
00435                         j++;
00436                 } else {
00437                         b--;
00438                 }
00439 
00440                 v_obs += bit_this ^ bit_prev;
00441 
00442                 bit_prev = bit_this;
00443         }
00444 
00445         v_obs += 1;
00446 
00447         numer = fabs(v_obs - (2.0 * numbits * pi * (1.0 - pi)));
00448         denom = 2.0 * sqrt(2.0 * numbits) * pi * (1.0 - pi);
00449 
00450         p_value = erfc(numer / denom);
00451 
00452         return (p_value);
00453 }
00454 
00455 /*
00456  * This is the block frequency test taken from the NIST SP 800-22 RNG
00457  * test suite.
00458  */
00459 static double
00460 blockfrequency(isc_mem_t *mctx, isc_uint16_t *values, size_t length) {
00461         isc_uint32_t i;
00462         isc_uint32_t numbits;
00463         isc_uint32_t mbits;
00464         isc_uint32_t mwords;
00465         isc_uint32_t numblocks;
00466         double *pi;
00467         isc_uint32_t cur_word;
00468         double chi_square;
00469         double p_value;
00470 
00471         numbits = length * 16;
00472         mbits = 32000;
00473         mwords = mbits / 16;
00474         numblocks = numbits / mbits;
00475 
00476         /* Debug message, not displayed when running via atf-run */
00477         printf("numblocks=%u\n", numblocks);
00478 
00479         /* Preconditions (section 2.2.7 in NIST SP 800-22) */
00480         ATF_REQUIRE(numbits >= 100);
00481         ATF_REQUIRE(mbits >= 20);
00482         ATF_REQUIRE((double) mbits > (0.01 * numbits));
00483         ATF_REQUIRE(numblocks < 100);
00484         ATF_REQUIRE(numbits >= (mbits * numblocks));
00485 
00486         pi = isc_mem_get(mctx, numblocks * sizeof(double));
00487         ATF_REQUIRE(pi != NULL);
00488 
00489         cur_word = 0;
00490         for (i = 0; i < numblocks; i++) {
00491                 isc_uint32_t j;
00492                 pi[i] = 0.0;
00493                 for (j = 0; j < mwords; j++) {
00494                         isc_uint32_t idx;
00495 
00496                         idx = i * mwords + j;
00497                         pi[i] += bitcounts_table[values[idx]];
00498                         cur_word++;
00499                 }
00500                 pi[i] /= mbits;
00501         }
00502 
00503         /* Compute chi_square */
00504         chi_square = 0.0;
00505         for (i = 0; i < numblocks; i++)
00506                 chi_square += (pi[i] - 0.5) * (pi[i] - 0.5);
00507 
00508         chi_square *= 4 * mbits;
00509 
00510         isc_mem_put(mctx, pi, numblocks * sizeof(double));
00511 
00512         /* Debug message, not displayed when running via atf-run */
00513         printf("chi_square=%f\n", chi_square);
00514 
00515         p_value = igamc(numblocks * 0.5, chi_square * 0.5);
00516 
00517         return (p_value);
00518 }
00519 
00520 /*
00521  * This is the binary matrix rank test taken from the NIST SP 800-22 RNG
00522  * test suite.
00523  */
00524 static double
00525 binarymatrixrank(isc_mem_t *mctx, isc_uint16_t *values, size_t length) {
00526         isc_uint32_t i;
00527         size_t matrix_m;
00528         size_t matrix_q;
00529         isc_uint32_t num_matrices;
00530         size_t numbits;
00531         isc_uint32_t fm_0;
00532         isc_uint32_t fm_1;
00533         isc_uint32_t fm_rest;
00534         double term1;
00535         double term2;
00536         double term3;
00537         double chi_square;
00538         double p_value;
00539 
00540         UNUSED(mctx);
00541 
00542         matrix_m = 32;
00543         matrix_q = 32;
00544         num_matrices = length / ((matrix_m * matrix_q) / 16);
00545         numbits = num_matrices * matrix_m * matrix_q;
00546 
00547         /* Preconditions (section 2.5.7 in NIST SP 800-22) */
00548         ATF_REQUIRE(matrix_m == 32);
00549         ATF_REQUIRE(matrix_q == 32);
00550         ATF_REQUIRE(numbits >= (38 * matrix_m * matrix_q));
00551 
00552         fm_0 = 0;
00553         fm_1 = 0;
00554         fm_rest = 0;
00555         for (i = 0; i < num_matrices; i++) {
00556                 /*
00557                  * Each isc_uint32_t supplies 32 bits, so a 32x32 bit matrix
00558                  * takes up isc_uint32_t array of size 32.
00559                  */
00560                 isc_uint32_t bits[32];
00561                 int j;
00562                 isc_uint32_t rank;
00563 
00564                 for (j = 0; j < 32; j++) {
00565                         size_t idx;
00566                         isc_uint32_t r1;
00567                         isc_uint32_t r2;
00568 
00569                         idx = i * ((matrix_m * matrix_q) / 16);
00570                         idx += j * 2;
00571 
00572                         r1 = values[idx];
00573                         r2 = values[idx + 1];
00574                         bits[j] = (r1 << 16) | r2;
00575                 }
00576 
00577                 rank = matrix_binaryrank(bits, matrix_m, matrix_q);
00578 
00579                 if (rank == matrix_m)
00580                         fm_0++;
00581                 else if (rank == (matrix_m - 1))
00582                         fm_1++;
00583                 else
00584                         fm_rest++;
00585         }
00586 
00587         /* Compute chi_square */
00588         term1 = ((fm_0 - (0.2888 * num_matrices)) *
00589                  (fm_0 - (0.2888 * num_matrices))) / (0.2888 * num_matrices);
00590         term2 = ((fm_1 - (0.5776 * num_matrices)) *
00591                  (fm_1 - (0.5776 * num_matrices))) / (0.5776 * num_matrices);
00592         term3 = ((fm_rest - (0.1336 * num_matrices)) *
00593                  (fm_rest - (0.1336 * num_matrices))) / (0.1336 * num_matrices);
00594 
00595         chi_square = term1 + term2 + term3;
00596 
00597         /* Debug message, not displayed when running via atf-run */
00598         printf("fm_0=%u, fm_1=%u, fm_rest=%u, chi_square=%f\n",
00599                fm_0, fm_1, fm_rest, chi_square);
00600 
00601         p_value = exp(-chi_square * 0.5);
00602 
00603         return (p_value);
00604 }
00605 
00606 ATF_TC(isc_rng_monobit);
00607 ATF_TC_HEAD(isc_rng_monobit, tc) {
00608         atf_tc_set_md_var(tc, "descr", "Monobit test for the RNG");
00609 }
00610 
00611 ATF_TC_BODY(isc_rng_monobit, tc) {
00612         UNUSED(tc);
00613 
00614         random_test(monobit);
00615 }
00616 
00617 ATF_TC(isc_rng_runs);
00618 ATF_TC_HEAD(isc_rng_runs, tc) {
00619         atf_tc_set_md_var(tc, "descr", "Runs test for the RNG");
00620 }
00621 
00622 ATF_TC_BODY(isc_rng_runs, tc) {
00623         UNUSED(tc);
00624 
00625         random_test(runs);
00626 }
00627 
00628 ATF_TC(isc_rng_blockfrequency);
00629 ATF_TC_HEAD(isc_rng_blockfrequency, tc) {
00630         atf_tc_set_md_var(tc, "descr", "Block frequency test for the RNG");
00631 }
00632 
00633 ATF_TC_BODY(isc_rng_blockfrequency, tc) {
00634         UNUSED(tc);
00635 
00636         random_test(blockfrequency);
00637 }
00638 
00639 ATF_TC(isc_rng_binarymatrixrank);
00640 ATF_TC_HEAD(isc_rng_binarymatrixrank, tc) {
00641         atf_tc_set_md_var(tc, "descr", "Binary matrix rank test for the RNG");
00642 }
00643 
00644 /*
00645  * This is the binary matrix rank test taken from the NIST SP 800-22 RNG
00646  * test suite.
00647  */
00648 ATF_TC_BODY(isc_rng_binarymatrixrank, tc) {
00649         UNUSED(tc);
00650 
00651         random_test(binarymatrixrank);
00652 }
00653 
00654 /*
00655  * Main
00656  */
00657 ATF_TP_ADD_TCS(tp) {
00658         ATF_TP_ADD_TC(tp, isc_rng_monobit);
00659         ATF_TP_ADD_TC(tp, isc_rng_runs);
00660         ATF_TP_ADD_TC(tp, isc_rng_blockfrequency);
00661         ATF_TP_ADD_TC(tp, isc_rng_binarymatrixrank);
00662 
00663         return (atf_no_error());
00664 }

Generated on Tue Apr 28 17:41:05 2015 by Doxygen 1.5.4 for BIND9 Internals 9.11.0pre-alpha