Home | History | Annotate | Download | only in n2rng
      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