1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 /* 22 * Copyright 2007 Sun Microsystems, Inc. All rights reserved. 23 * Use is subject to license terms. 24 */ 25 26 #pragma ident "%Z%%M% %I% %E% SMI" 27 28 29 #include <sys/ddi.h> 30 #include <sys/errno.h> 31 #include <sys/types.h> 32 #include <sys/n2rng.h> 33 #include <sys/int_types.h> 34 35 36 /* 37 * This whole file is really doing floating point type stuff, and 38 * would be quite simple in user space. But since we are in the 39 * kernel, (a) we can't use floating point, and (b) we don't have a 40 * math library. 41 */ 42 43 /* used inside msb */ 44 #define MSBSTEP(word, shift, counter) \ 45 if (word & (~0ULL << shift)) { \ 46 word >>= shift; \ 47 counter += shift; \ 48 } 49 50 /* 51 * returns the position of the MSB of x. The 1 bit is position 0. An 52 * all zero arg returns -1. 53 */ 54 static int 55 msb(uint64_t x) 56 { 57 int bit; 58 59 if (x == 0) { 60 return (-1); 61 } 62 63 bit = 0; 64 MSBSTEP(x, 32, bit); 65 MSBSTEP(x, 16, bit); 66 MSBSTEP(x, 8, bit); 67 MSBSTEP(x, 4, bit); 68 MSBSTEP(x, 2, bit); 69 MSBSTEP(x, 1, bit); 70 71 return (bit); 72 } 73 74 /* 75 * lg2 computes 2^(LOG_VAL_SCALE) * log2(x/2^LOG_ARG_SCALE), where ^ 76 * is exponentiation. 77 * 78 * The following conditions must be satisfied: LOG_VAL_SCALE <= 62, 79 * LOG_VAL_SCALE + log2(maxarg) < 64, LOG_VAL_SCALE >= 0, 80 * LOG_ARG_SCALE <= 63. Recommended LOG_VAL_SCALE is 57, which is the 81 * largest value such that overflow is impossible. 82 */ 83 static int64_t 84 lg2(uint64_t x) 85 { 86 /* 87 * logtable[i-1] == round(2^63 * log2(2^i/(2^i - 1))), where ^ 88 * is exponentiation. 89 */ 90 static const uint64_t logtable[] = { 91 9223372036854775808ULL, 3828045265094622256ULL, 92 1776837224931603046ULL, 858782676832593460ULL, 93 422464469962470743ULL, 209555718266071751ULL, 94 104365343613858422ULL, 52080352580344565ULL, 95 26014696649359209ULL, 13000990870918027ULL, 96 6498907625079429ULL, 3249057053828501ULL, 97 1624429361456373ULL, 812189892390238ULL, 98 406088749488886ULL, 203042825615163ULL, 99 101521025531171ULL, 50760415947221ULL, 100 25380183769112ULL, 12690085833443ULL, 101 6345041403945ULL, 3172520323778ULL, 102 1586260067341ULL, 793130010033ULL, 103 396564999107ULL, 198282498076ULL, 104 99141248669ULL, 49570624242ULL, 105 24785312098ULL, 12392656043ULL, 6196328020ULL, 3098164010ULL, 106 1549082005ULL, 774541002ULL, 387270501ULL, 193635251ULL, 107 96817625ULL, 48408813ULL, 24204406ULL, 12102203ULL, 6051102ULL, 108 3025551ULL, 1512775ULL, 756388ULL, 378194ULL, 189097ULL, 109 94548ULL, 47274ULL, 23637ULL, 11819ULL, 5909ULL, 2955ULL, 110 1477ULL, 739ULL, 369ULL, 185ULL, 92ULL, 46ULL, 23ULL, 111 12ULL, 6ULL, 3ULL, 1ULL 112 }; 113 114 uint64_t xx; 115 uint64_t logx; 116 uint64_t tmp; 117 int i; 118 119 if (x == 0) { 120 return (-INT64_MAX - 1); 121 } 122 123 /* 124 * Invariant: log2(xx) + logx == log2(x). This is true at the after 125 * the normalization. At each adjustment step we multiply xx by 126 * (2^i-1)/2^i, which effectively decreases log2(xx) by 127 * log2(2^i/(2^i-1)), and a the same time, we add table[i], which 128 * equals log2(2^i/(2^i-1)), to logx. By induction the invariant is 129 * true at the end. At the end xx==1, so log2(xx)==0, and thus 130 * logx=log2(x); 131 */ 132 /* Normalize */ 133 i = msb(x); /* use i in computing preshift */ 134 if (i - LOG_ARG_SCALE > 0) { 135 xx = x >> (i - LOG_ARG_SCALE); 136 } else { 137 xx = x << (LOG_ARG_SCALE - i); 138 } 139 logx = (int64_t)(i - LOG_ARG_SCALE) << LOG_VAL_SCALE; 140 141 for (i = 1; i <= LOG_ARG_SCALE; i++) { 142 /* 1ULL << (i-1) is rounding */ 143 while ((tmp = xx - ((xx + (1ULL << (i-1))) >> i)) >= 144 1ULL << LOG_ARG_SCALE) { 145 xx = tmp; 146 /* 1ULL << (63 - LOG_VAL_SCALE -1) is rounding */ 147 logx += (logtable[i-1] + 148 (1ULL << (63 - LOG_VAL_SCALE - 1))) >> 149 (63 - LOG_VAL_SCALE); 150 } 151 } 152 153 return (logx); 154 } 155 156 157 158 /* 159 * The EXCHANGE macro swaps entries j & k if necessary so that 160 * data[j] <= data[k]. 161 * 162 * If OBLIVIOUS is defined, no branches are used. This would allow 163 * this algorithm to be used by the CPU manufacturing people who run 164 * on a tester that requires the exact same instruction address stream 165 * on every test. (It's a bit slower with OBLIVIOUS defined.) 166 */ 167 #ifdef OBLIVIOUS 168 #define EXCHANGE(j, k) \ 169 { \ 170 uint64_t tmp, mask; \ 171 mask = (uint64_t)(((int64_t)(data[k] - data[j])) >> 63); \ 172 tmp = data[j] + data[k]; \ 173 data[j] = data[k] & mask | data[j] & ~mask; \ 174 data[k] = tmp - data[j]; \ 175 } 176 #else 177 #define EXCHANGE(j, k) \ 178 { \ 179 uint64_t tmp; \ 180 if (data[j] > data[k]) { \ 181 tmp = data[j]; \ 182 data[j] = data[k]; \ 183 data[k] = tmp; \ 184 } \ 185 } 186 #endif 187 188 189 190 /* 191 * This is a Batcher sort from Knuth v. 3. There is no flow control 192 * that depends on the values being sorted, except in the EXCHANGE 193 * step, but that can be made oblivious to the data values, too, by 194 * setting OBLIVIOUS. So this code could be using in chip testers 195 * that require fixed flow through a test. 196 * 197 * This is presently hard-coded for sorting uint64_t values. 198 */ 199 void 200 n2rng_sort(uint64_t *data, int log2_size) 201 { 202 int p, q, d, r, i; 203 204 for (p = 1 << (log2_size - 1); p > 0; p >>= 1) { 205 d = p; 206 r = 0; 207 for (q = 1 << (log2_size - 1); q >= p; q >>= 1) { 208 for (i = 0; i + d < (1 << log2_size); i++) { 209 if ((i & p) == r) { 210 EXCHANGE(i, i+d); 211 } 212 } 213 d = q - p; 214 r = p; 215 } 216 } 217 } 218 219 220 /* 221 * Computes several measures of entropy per word: Renyi H0 (log2 of 222 * number of distinct symbols), Renyi H1 (Shannon), 223 * Renyi H2 (-log2 of sum(P_i^2)), and 224 * Renyi H-infinity (min). The results are coded as H * 225 * 2^LOG_VAL_SCALE). The samples array is modified by sorting in 226 * place. 227 * 228 * None if this is really valid, since it requres that the block 229 * length be at least as long as the largest non-approximately-zero 230 * coefficient in the autocorrelation function, and that the number 231 * of samples be much larger than 2^longest_block_length_in_bits. 232 * But we hope that bigger is better, even when it is invalid. 233 */ 234 void 235 n2rng_renyi_entropy(uint64_t *samples, int lg2samples, n2rng_osc_perf_t *entp) 236 { 237 size_t i; 238 uint64_t cv = samples[0]; /* current value */ 239 size_t count = 1; 240 size_t numdistinct = 0; 241 size_t largestcount = 0; 242 uint64_t shannonsum = 0; 243 uint64_t sqsum = 0; 244 245 n2rng_sort(samples, lg2samples); 246 247 for (i = 1; i < (1 << lg2samples); i++) { 248 if (samples[i] != cv) { 249 numdistinct++; 250 if (count > largestcount) { 251 largestcount = count; 252 } 253 #ifdef COMPUTE_SHANNON_ENTROPY 254 shannonsum -= (count * (lg2(count) + 255 ((int64_t)(LOG_ARG_SCALE - lg2samples) << 256 LOG_VAL_SCALE))) >> lg2samples; 257 #endif /* COMPUTE_SHANNON_ENTROPY */ 258 sqsum += count * count; 259 count = 1; 260 cv = samples[i]; 261 } else { 262 count++; 263 } 264 } 265 /* process last block */ 266 numdistinct++; 267 if (count > largestcount) { 268 largestcount = count; 269 } 270 #ifdef COMPUTE_SHANNON_ENTROPY 271 shannonsum -= (count * (lg2(count) + 272 ((int64_t)(LOG_ARG_SCALE - lg2samples) << LOG_VAL_SCALE))) >> 273 lg2samples; 274 #endif /* COMPUTE_SHANNON_ENTROPY */ 275 sqsum += count * count; 276 277 entp->numvals = numdistinct; 278 /* H1 is shannon entropy: -sum(p_i * log2(p_i)) */ 279 entp->H1 = shannonsum / 64; 280 /* H2 is -log2(sum p_i^2) */ 281 entp->H2 = -(lg2(sqsum) + 282 ((int64_t)(LOG_ARG_SCALE - 2 * lg2samples) << LOG_VAL_SCALE)) / 64; 283 /* Hinf = -log2(highest_probability) */ 284 entp->Hinf = -(lg2(largestcount) + 285 ((int64_t)(LOG_ARG_SCALE - lg2samples) << LOG_VAL_SCALE)) / 64; 286 } 287