mirror of
https://github.com/jart/cosmopolitan.git
synced 2025-07-04 10:18:31 +00:00
Fold LIBC_RAND into LIBC_STDIO/TINYMATH/INTRIN
This commit is contained in:
Justine Tunney
parent
05b8f82371
commit
8a0a2c0c36
183 changed files with 149 additions and 322 deletions
54
libc/tinymath/measureentropy.c
Normal file
54
libc/tinymath/measureentropy.c
Normal file
|
|
@ -0,0 +1,54 @@
|
|||
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
|
||||
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
|
||||
╞══════════════════════════════════════════════════════════════════════════════╡
|
||||
│ Copyright 2021 Justine Alexandra Roberts Tunney │
|
||||
│ │
|
||||
│ Permission to use, copy, modify, and/or distribute this software for │
|
||||
│ any purpose with or without fee is hereby granted, provided that the │
|
||||
│ above copyright notice and this permission notice appear in all copies. │
|
||||
│ │
|
||||
│ THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL │
|
||||
│ WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED │
|
||||
│ WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE │
|
||||
│ AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL │
|
||||
│ DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR │
|
||||
│ PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER │
|
||||
│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │
|
||||
│ PERFORMANCE OF THIS SOFTWARE. │
|
||||
╚─────────────────────────────────────────────────────────────────────────────*/
|
||||
#include "libc/math.h"
|
||||
#include "libc/stdio/rand.h"
|
||||
#include "libc/str/str.h"
|
||||
|
||||
/**
|
||||
* Returns Shannon entropy of array.
|
||||
*
|
||||
* This gives you an idea of the density of information. Cryptographic
|
||||
* random should be in the ballpark of 7.9 whereas plaintext will be
|
||||
* more like 4.5.
|
||||
*
|
||||
* @param p is treated as binary octets
|
||||
* @param n should be at least 1000
|
||||
* @return number between 0 and 8
|
||||
*/
|
||||
double MeasureEntropy(const char *p, size_t n) {
|
||||
size_t i;
|
||||
double e, x;
|
||||
long h[256];
|
||||
e = 0;
|
||||
if (n) {
|
||||
bzero(h, sizeof(h));
|
||||
for (i = 0; i < n; ++i) {
|
||||
++h[p[i] & 255];
|
||||
}
|
||||
for (i = 0; i < 256; i++) {
|
||||
if (h[i]) {
|
||||
x = h[i];
|
||||
x /= n;
|
||||
e += x * log(x);
|
||||
}
|
||||
}
|
||||
e = -(e / M_LN2);
|
||||
}
|
||||
return e;
|
||||
}
|
||||
142
libc/tinymath/poz.c
Normal file
142
libc/tinymath/poz.c
Normal file
|
|
@ -0,0 +1,142 @@
|
|||
/*-*- mode:c;indent-tabs-mode:t;c-basic-offset:4;tab-width:4;coding:utf-8 -*-│
|
||||
│vi: set et ft=c ts=4 sts=4 sw=4 fenc=utf-8 :vi│
|
||||
╚─────────────────────────────────────────────────────────────────────────────*/
|
||||
/* clang-format off */
|
||||
/*
|
||||
|
||||
Compute probability of measured Chi Square value.
|
||||
|
||||
This code was developed by Gary Perlman of the Wang
|
||||
Institute (full citation below) and has been minimally
|
||||
modified for use in this program.
|
||||
|
||||
*/
|
||||
|
||||
#include "libc/math.h"
|
||||
|
||||
/*HEADER
|
||||
Module: z.c
|
||||
Purpose: compute approximations to normal z distribution probabilities
|
||||
Programmer: Gary Perlman
|
||||
Organization: Wang Institute, Tyngsboro, MA 01879
|
||||
Copyright: none
|
||||
Tabstops: 4
|
||||
*/
|
||||
|
||||
#define Z_MAX 6.0 /* maximum meaningful z value */
|
||||
|
||||
/*FUNCTION poz: probability of normal z value */
|
||||
/*ALGORITHM
|
||||
Adapted from a polynomial approximation in:
|
||||
Ibbetson D, Algorithm 209
|
||||
Collected Algorithms of the CACM 1963 p. 616
|
||||
Note:
|
||||
This routine has six digit accuracy, so it is only useful for absolute
|
||||
z values < 6. For z values >= to 6.0, poz() returns 0.0.
|
||||
*/
|
||||
static double /*VAR returns cumulative probability from -oo to z */
|
||||
poz(const double z) /*VAR normal z value */
|
||||
{
|
||||
double y, x, w;
|
||||
|
||||
if (z == 0.0) {
|
||||
x = 0.0;
|
||||
} else {
|
||||
y = 0.5 * fabs(z);
|
||||
if (y >= (Z_MAX * 0.5)) {
|
||||
x = 1.0;
|
||||
} else if (y < 1.0) {
|
||||
w = y * y;
|
||||
x = ((((((((0.000124818987 * w
|
||||
-0.001075204047) * w +0.005198775019) * w
|
||||
-0.019198292004) * w +0.059054035642) * w
|
||||
-0.151968751364) * w +0.319152932694) * w
|
||||
-0.531923007300) * w +0.797884560593) * y * 2.0;
|
||||
} else {
|
||||
y -= 2.0;
|
||||
x = (((((((((((((-0.000045255659 * y
|
||||
+0.000152529290) * y -0.000019538132) * y
|
||||
-0.000676904986) * y +0.001390604284) * y
|
||||
-0.000794620820) * y -0.002034254874) * y
|
||||
+0.006549791214) * y -0.010557625006) * y
|
||||
+0.011630447319) * y -0.009279453341) * y
|
||||
+0.005353579108) * y -0.002141268741) * y
|
||||
+0.000535310849) * y +0.999936657524;
|
||||
}
|
||||
}
|
||||
return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
|
||||
}
|
||||
|
||||
/*
|
||||
Module: chisq.c
|
||||
Purpose: compute approximations to chisquare distribution probabilities
|
||||
Contents: pochisq()
|
||||
Uses: poz() in z.c (Algorithm 209)
|
||||
Programmer: Gary Perlman
|
||||
Organization: Wang Institute, Tyngsboro, MA 01879
|
||||
Copyright: none
|
||||
Tabstops: 4
|
||||
*/
|
||||
|
||||
#define LOG_SQRT_PI 0.5723649429247000870717135 /* log (sqrt (pi)) */
|
||||
#define I_SQRT_PI 0.5641895835477562869480795 /* 1 / sqrt (pi) */
|
||||
#define BIGX 20.0 /* max value to represent exp (x) */
|
||||
#define ex(x) (((x) < -BIGX) ? 0.0 : exp(x))
|
||||
|
||||
/*FUNCTION pochisq: probability of chi square value */
|
||||
/*ALGORITHM Compute probability of chi square value.
|
||||
Adapted from:
|
||||
Hill, I. D. and Pike, M. C. Algorithm 299
|
||||
Collected Algorithms for the CACM 1967 p. 243
|
||||
Updated for rounding errors based on remark in
|
||||
ACM TOMS June 1985, page 185
|
||||
*/
|
||||
|
||||
double pochisq(
|
||||
const double ax, /* obtained chi-square value */
|
||||
const int df /* degrees of freedom */
|
||||
)
|
||||
{
|
||||
double x = ax;
|
||||
double a, y, s;
|
||||
double e, c, z;
|
||||
int even; /* true if df is an even number */
|
||||
|
||||
y = 0.0; /* XXX: blind modification due to GCC error */
|
||||
|
||||
if (x <= 0.0 || df < 1) {
|
||||
return 1.0;
|
||||
}
|
||||
|
||||
a = 0.5 * x;
|
||||
even = (2 * (df / 2)) == df;
|
||||
if (df > 1) {
|
||||
y = ex(-a);
|
||||
}
|
||||
s = (even ? y : (2.0 * poz(-sqrt(x))));
|
||||
if (df > 2) {
|
||||
x = 0.5 * (df - 1.0);
|
||||
z = (even ? 1.0 : 0.5);
|
||||
if (a > BIGX) {
|
||||
e = (even ? 0.0 : LOG_SQRT_PI);
|
||||
c = log(a);
|
||||
while (z <= x) {
|
||||
e = log(z) + e;
|
||||
s += ex(c * z - a - e);
|
||||
z += 1.0;
|
||||
}
|
||||
return (s);
|
||||
} else {
|
||||
e = (even ? 1.0 : (I_SQRT_PI / sqrt(a)));
|
||||
c = 0.0;
|
||||
while (z <= x) {
|
||||
e = e * (a / z);
|
||||
c = c + e;
|
||||
z += 1.0;
|
||||
}
|
||||
return (c * y + s);
|
||||
}
|
||||
} else {
|
||||
return s;
|
||||
}
|
||||
}
|
||||
184
libc/tinymath/randtest.c
Normal file
184
libc/tinymath/randtest.c
Normal file
|
|
@ -0,0 +1,184 @@
|
|||
/* clang-format off */
|
||||
/*
|
||||
|
||||
Apply various randomness tests to a stream of bytes
|
||||
|
||||
by John Walker -- September 1996
|
||||
http://www.fourmilab.ch/
|
||||
|
||||
*/
|
||||
|
||||
#include "libc/math.h"
|
||||
|
||||
#define FALSE 0
|
||||
#define TRUE 1
|
||||
|
||||
#define log2of10 3.32192809488736234787
|
||||
|
||||
static int binary = FALSE; /* Treat input as a bitstream */
|
||||
|
||||
static long ccount[256], /* Bins to count occurrences of values */
|
||||
totalc = 0; /* Total bytes counted */
|
||||
static double prob[256]; /* Probabilities per bin for entropy */
|
||||
|
||||
/* RT_LOG2 -- Calculate log to the base 2 */
|
||||
|
||||
static double rt_log2(double x)
|
||||
{
|
||||
return log2of10 * log10(x);
|
||||
}
|
||||
|
||||
#define MONTEN 6 /* Bytes used as Monte Carlo
|
||||
co-ordinates. This should be no more
|
||||
bits than the mantissa of your
|
||||
"double" floating point type. */
|
||||
|
||||
static int mp, sccfirst;
|
||||
static unsigned int monte[MONTEN];
|
||||
static long inmont, mcount_;
|
||||
static double cexp_, incirc, montex, montey, montepi,
|
||||
scc, sccun, sccu0, scclast, scct1, scct2, scct3,
|
||||
ent, chisq, datasum;
|
||||
|
||||
/* RT_INIT -- Initialise random test counters. */
|
||||
|
||||
void rt_init(int binmode)
|
||||
{
|
||||
int i;
|
||||
|
||||
binary = binmode; /* Set binary / byte mode */
|
||||
|
||||
/* Initialise for calculations */
|
||||
|
||||
ent = 0.0; /* Clear entropy accumulator */
|
||||
chisq = 0.0; /* Clear Chi-Square */
|
||||
datasum = 0.0; /* Clear sum of bytes for arithmetic mean */
|
||||
|
||||
mp = 0; /* Reset Monte Carlo accumulator pointer */
|
||||
mcount_ = 0; /* Clear Monte Carlo tries */
|
||||
inmont = 0; /* Clear Monte Carlo inside count */
|
||||
incirc = 65535.0 * 65535.0;/* In-circle distance for Monte Carlo */
|
||||
|
||||
sccfirst = TRUE; /* Mark first time for serial correlation */
|
||||
scct1 = scct2 = scct3 = 0.0; /* Clear serial correlation terms */
|
||||
|
||||
incirc = pow(pow(256.0, (double) (MONTEN / 2)) - 1, 2.0);
|
||||
|
||||
for (i = 0; i < 256; i++) {
|
||||
ccount[i] = 0;
|
||||
}
|
||||
totalc = 0;
|
||||
}
|
||||
|
||||
/* RT_ADD -- Add one or more bytes to accumulation. */
|
||||
|
||||
void rt_add(void *buf, int bufl)
|
||||
{
|
||||
unsigned char *bp = buf;
|
||||
int oc, c, bean;
|
||||
|
||||
while (bean = 0, (bufl-- > 0)) {
|
||||
oc = *bp++;
|
||||
|
||||
do {
|
||||
if (binary) {
|
||||
c = !!(oc & 0x80);
|
||||
} else {
|
||||
c = oc;
|
||||
}
|
||||
ccount[c]++; /* Update counter for this bin */
|
||||
totalc++;
|
||||
|
||||
/* Update inside / outside circle counts for Monte Carlo
|
||||
computation of PI */
|
||||
|
||||
if (bean == 0) {
|
||||
monte[mp++] = oc; /* Save character for Monte Carlo */
|
||||
if (mp >= MONTEN) { /* Calculate every MONTEN character */
|
||||
int mj;
|
||||
|
||||
mp = 0;
|
||||
mcount_++;
|
||||
montex = montey = 0;
|
||||
for (mj = 0; mj < MONTEN / 2; mj++) {
|
||||
montex = (montex * 256.0) + monte[mj];
|
||||
montey = (montey * 256.0) + monte[(MONTEN / 2) + mj];
|
||||
}
|
||||
if ((montex * montex + montey * montey) <= incirc) {
|
||||
inmont++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Update calculation of serial correlation coefficient */
|
||||
|
||||
sccun = c;
|
||||
if (sccfirst) {
|
||||
sccfirst = FALSE;
|
||||
scclast = 0;
|
||||
sccu0 = sccun;
|
||||
} else {
|
||||
scct1 = scct1 + scclast * sccun;
|
||||
}
|
||||
scct2 = scct2 + sccun;
|
||||
scct3 = scct3 + (sccun * sccun);
|
||||
scclast = sccun;
|
||||
oc <<= 1;
|
||||
} while (binary && (++bean < 8));
|
||||
}
|
||||
}
|
||||
|
||||
/* RT_END -- Complete calculation and return results. */
|
||||
|
||||
void rt_end(double *r_ent, double *r_chisq, double *r_mean,
|
||||
double *r_montepicalc, double *r_scc)
|
||||
{
|
||||
int i;
|
||||
|
||||
/* Complete calculation of serial correlation coefficient */
|
||||
|
||||
scct1 = scct1 + scclast * sccu0;
|
||||
scct2 = scct2 * scct2;
|
||||
scc = totalc * scct3 - scct2;
|
||||
if (scc == 0.0) {
|
||||
scc = -100000;
|
||||
} else {
|
||||
scc = (totalc * scct1 - scct2) / scc;
|
||||
}
|
||||
|
||||
/* Scan bins and calculate probability for each bin and
|
||||
Chi-Square distribution. The probability will be reused
|
||||
in the entropy calculation below. While we're at it,
|
||||
we sum of all the data which will be used to compute the
|
||||
mean. */
|
||||
|
||||
cexp_ = totalc / (binary ? 2.0 : 256.0); /* Expected count per bin */
|
||||
for (i = 0; i < (binary ? 2 : 256); i++) {
|
||||
double a = ccount[i] - cexp_;
|
||||
|
||||
prob[i] = ((double) ccount[i]) / totalc;
|
||||
chisq += (a * a) / cexp_;
|
||||
datasum += ((double) i) * ccount[i];
|
||||
}
|
||||
|
||||
/* Calculate entropy */
|
||||
|
||||
for (i = 0; i < (binary ? 2 : 256); i++) {
|
||||
if (prob[i] > 0.0) {
|
||||
ent += prob[i] * rt_log2(1 / prob[i]);
|
||||
}
|
||||
}
|
||||
|
||||
/* Calculate Monte Carlo value for PI from percentage of hits
|
||||
within the circle */
|
||||
|
||||
montepi = 4.0 * (((double) inmont) / mcount_);
|
||||
|
||||
/* Return results through arguments */
|
||||
|
||||
*r_ent = ent;
|
||||
*r_chisq = chisq;
|
||||
*r_mean = datasum / totalc;
|
||||
*r_montepicalc = montepi;
|
||||
*r_scc = scc;
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue