/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License (the "License").
* You may not use this file except in compliance with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* or http://www.opensolaris.org/os/licensing.
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2007 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
#pragma ident "%Z%%M% %I% %E% SMI"
#include <sys/ddi.h>
#include <sys/errno.h>
#include <sys/types.h>
#include <sys/n2rng.h>
#include <sys/int_types.h>
/*
* This whole file is really doing floating point type stuff, and
* would be quite simple in user space. But since we are in the
* kernel, (a) we can't use floating point, and (b) we don't have a
* math library.
*/
/* used inside msb */
#define MSBSTEP(word, shift, counter) \
if (word & (~0ULL << shift)) { \
word >>= shift; \
counter += shift; \
}
/*
* returns the position of the MSB of x. The 1 bit is position 0. An
* all zero arg returns -1.
*/
static int
msb(uint64_t x)
{
int bit;
if (x == 0) {
return (-1);
}
bit = 0;
MSBSTEP(x, 32, bit);
MSBSTEP(x, 16, bit);
MSBSTEP(x, 8, bit);
MSBSTEP(x, 4, bit);
MSBSTEP(x, 2, bit);
MSBSTEP(x, 1, bit);
return (bit);
}
/*
* lg2 computes 2^(LOG_VAL_SCALE) * log2(x/2^LOG_ARG_SCALE), where ^
* is exponentiation.
*
* The following conditions must be satisfied: LOG_VAL_SCALE <= 62,
* LOG_VAL_SCALE + log2(maxarg) < 64, LOG_VAL_SCALE >= 0,
* LOG_ARG_SCALE <= 63. Recommended LOG_VAL_SCALE is 57, which is the
* largest value such that overflow is impossible.
*/
static int64_t
lg2(uint64_t x)
{
/*
* logtable[i-1] == round(2^63 * log2(2^i/(2^i - 1))), where ^
* is exponentiation.
*/
static const uint64_t logtable[] = {
9223372036854775808ULL, 3828045265094622256ULL,
1776837224931603046ULL, 858782676832593460ULL,
422464469962470743ULL, 209555718266071751ULL,
104365343613858422ULL, 52080352580344565ULL,
26014696649359209ULL, 13000990870918027ULL,
6498907625079429ULL, 3249057053828501ULL,
1624429361456373ULL, 812189892390238ULL,
406088749488886ULL, 203042825615163ULL,
101521025531171ULL, 50760415947221ULL,
25380183769112ULL, 12690085833443ULL,
6345041403945ULL, 3172520323778ULL,
1586260067341ULL, 793130010033ULL,
396564999107ULL, 198282498076ULL,
99141248669ULL, 49570624242ULL,
24785312098ULL, 12392656043ULL, 6196328020ULL, 3098164010ULL,
1549082005ULL, 774541002ULL, 387270501ULL, 193635251ULL,
96817625ULL, 48408813ULL, 24204406ULL, 12102203ULL, 6051102ULL,
3025551ULL, 1512775ULL, 756388ULL, 378194ULL, 189097ULL,
94548ULL, 47274ULL, 23637ULL, 11819ULL, 5909ULL, 2955ULL,
1477ULL, 739ULL, 369ULL, 185ULL, 92ULL, 46ULL, 23ULL,
12ULL, 6ULL, 3ULL, 1ULL
};
uint64_t xx;
uint64_t logx;
uint64_t tmp;
int i;
if (x == 0) {
return (-INT64_MAX - 1);
}
/*
* Invariant: log2(xx) + logx == log2(x). This is true at the after
* the normalization. At each adjustment step we multiply xx by
* (2^i-1)/2^i, which effectively decreases log2(xx) by
* log2(2^i/(2^i-1)), and a the same time, we add table[i], which
* equals log2(2^i/(2^i-1)), to logx. By induction the invariant is
* true at the end. At the end xx==1, so log2(xx)==0, and thus
* logx=log2(x);
*/
/* Normalize */
i = msb(x); /* use i in computing preshift */
if (i - LOG_ARG_SCALE > 0) {
xx = x >> (i - LOG_ARG_SCALE);
} else {
xx = x << (LOG_ARG_SCALE - i);
}
logx = (int64_t)(i - LOG_ARG_SCALE) << LOG_VAL_SCALE;
for (i = 1; i <= LOG_ARG_SCALE; i++) {
/* 1ULL << (i-1) is rounding */
while ((tmp = xx - ((xx + (1ULL << (i-1))) >> i)) >=
1ULL << LOG_ARG_SCALE) {
xx = tmp;
/* 1ULL << (63 - LOG_VAL_SCALE -1) is rounding */
logx += (logtable[i-1] +
(1ULL << (63 - LOG_VAL_SCALE - 1))) >>
(63 - LOG_VAL_SCALE);
}
}
return (logx);
}
/*
* The EXCHANGE macro swaps entries j & k if necessary so that
* data[j] <= data[k].
*
* If OBLIVIOUS is defined, no branches are used. This would allow
* this algorithm to be used by the CPU manufacturing people who run
* on a tester that requires the exact same instruction address stream
* on every test. (It's a bit slower with OBLIVIOUS defined.)
*/
#ifdef OBLIVIOUS
#define EXCHANGE(j, k) \
{ \
uint64_t tmp, mask; \
mask = (uint64_t)(((int64_t)(data[k] - data[j])) >> 63); \
tmp = data[j] + data[k]; \
data[j] = data[k] & mask | data[j] & ~mask; \
data[k] = tmp - data[j]; \
}
#else
#define EXCHANGE(j, k) \
{ \
uint64_t tmp; \
if (data[j] > data[k]) { \
tmp = data[j]; \
data[j] = data[k]; \
data[k] = tmp; \
} \
}
#endif
/*
* This is a Batcher sort from Knuth v. 3. There is no flow control
* that depends on the values being sorted, except in the EXCHANGE
* step, but that can be made oblivious to the data values, too, by
* setting OBLIVIOUS. So this code could be using in chip testers
* that require fixed flow through a test.
*
* This is presently hard-coded for sorting uint64_t values.
*/
void
n2rng_sort(uint64_t *data, int log2_size)
{
int p, q, d, r, i;
for (p = 1 << (log2_size - 1); p > 0; p >>= 1) {
d = p;
r = 0;
for (q = 1 << (log2_size - 1); q >= p; q >>= 1) {
for (i = 0; i + d < (1 << log2_size); i++) {
if ((i & p) == r) {
EXCHANGE(i, i+d);
}
}
d = q - p;
r = p;
}
}
}
/*
* Computes several measures of entropy per word: Renyi H0 (log2 of
* number of distinct symbols), Renyi H1 (Shannon),
* Renyi H2 (-log2 of sum(P_i^2)), and
* Renyi H-infinity (min). The results are coded as H *
* 2^LOG_VAL_SCALE). The samples array is modified by sorting in
* place.
*
* None if this is really valid, since it requres that the block
* length be at least as long as the largest non-approximately-zero
* coefficient in the autocorrelation function, and that the number
* of samples be much larger than 2^longest_block_length_in_bits.
* But we hope that bigger is better, even when it is invalid.
*/
void
n2rng_renyi_entropy(uint64_t *samples, int lg2samples, n2rng_osc_perf_t *entp)
{
size_t i;
uint64_t cv = samples[0]; /* current value */
size_t count = 1;
size_t numdistinct = 0;
size_t largestcount = 0;
uint64_t shannonsum = 0;
uint64_t sqsum = 0;
n2rng_sort(samples, lg2samples);
for (i = 1; i < (1 << lg2samples); i++) {
if (samples[i] != cv) {
numdistinct++;
if (count > largestcount) {
largestcount = count;
}
#ifdef COMPUTE_SHANNON_ENTROPY
shannonsum -= (count * (lg2(count) +
((int64_t)(LOG_ARG_SCALE - lg2samples) <<
LOG_VAL_SCALE))) >> lg2samples;
#endif /* COMPUTE_SHANNON_ENTROPY */
sqsum += count * count;
count = 1;
cv = samples[i];
} else {
count++;
}
}
/* process last block */
numdistinct++;
if (count > largestcount) {
largestcount = count;
}
#ifdef COMPUTE_SHANNON_ENTROPY
shannonsum -= (count * (lg2(count) +
((int64_t)(LOG_ARG_SCALE - lg2samples) << LOG_VAL_SCALE))) >>
lg2samples;
#endif /* COMPUTE_SHANNON_ENTROPY */
sqsum += count * count;
entp->numvals = numdistinct;
/* H1 is shannon entropy: -sum(p_i * log2(p_i)) */
entp->H1 = shannonsum / 64;
/* H2 is -log2(sum p_i^2) */
entp->H2 = -(lg2(sqsum) +
((int64_t)(LOG_ARG_SCALE - 2 * lg2samples) << LOG_VAL_SCALE)) / 64;
/* Hinf = -log2(highest_probability) */
entp->Hinf = -(lg2(largestcount) +
((int64_t)(LOG_ARG_SCALE - lg2samples) << LOG_VAL_SCALE)) / 64;
}