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
msb(uint64_t x)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
lg2(uint64_t x)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
n2rng_sort(uint64_t * data,int log2_size)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
n2rng_renyi_entropy(uint64_t * samples,int lg2samples,n2rng_osc_perf_t * entp)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