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)  \
45if (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 */
54static int
55msb(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 */
83static int64_t
84lg2(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 */
199void
200n2rng_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 */
234void
235n2rng_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