Entropy Visualizer — C source

/*
 * entropy-visualizer — byte-level randomness analysis: Shannon entropy, a χ²
 *                      uniformity test, serial correlation and a Monte Carlo π
 *                      estimate.
 *
 * Language: C (C11, standard library only — <math.h> covers everything here)
 * Source:   CosmoDev polyglot showcase port of the Entropy Visualizer tool,
 *           ported from src/lib/entropy-visualizer.ts (the canonical TypeScript
 *           implementation).
 * License:  display source — part of CosmoDev's polyglot tool pages.
 *
 * Everything is deterministic: the same bytes always produce the same numbers.
 * The statistics are the classic `ent`-style battery — high Shannon entropy
 * alone does not prove randomness, which is why the χ² p-value, the serial
 * correlation and the π estimate are reported alongside it. Compressed data
 * scores ~7.99 bits/byte yet fails χ² badly.
 *
 * Build: cc -std=c11 entropy-visualizer.c -lm
 */

#include <math.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* --------------------------------------------------------------- constants --- */

/* Bytes per heatmap block. */
#define BLOCK_SIZE 16

/* Degrees of freedom for the byte-frequency χ² test (256 bins - 1 constraint). */
#define CHI_SQUARED_DF 255

typedef enum {
    VERDICT_EXCELLENT,
    VERDICT_GOOD,
    VERDICT_SUSPICIOUS,
    VERDICT_LOW
} entropy_verdict;

const char *verdict_name(entropy_verdict v)
{
    switch (v) {
        case VERDICT_EXCELLENT:  return "excellent";
        case VERDICT_GOOD:       return "good";
        case VERDICT_SUSPICIOUS: return "suspicious";
        default:                 return "low";
    }
}

typedef struct {
    /* Shannon entropy in bits per byte (0 = one repeating byte, 8 = uniform). */
    double shannon_entropy;
    /* χ² statistic against the uniform 256-bin expectation. */
    double chi_squared;
    /* Upper-tail p-value for χ² with 255 df (1 = perfectly plausible). */
    double chi_squared_p_value;
    /* Circular serial correlation between consecutive bytes (-1 … +1). */
    double serial_correlation;
    /* Monte Carlo π estimate from consecutive byte pairs (≈3.14159 if random). */
    double monte_carlo_pi;
    /* Occurrence count per byte value 0-255. */
    long byte_frequencies[256];
    /* Shannon entropy of each BLOCK_SIZE-byte block, for the heatmap. */
    double *block_entropies;
    size_t  block_count;
    entropy_verdict verdict;
} entropy_analysis;

/* -------------------------------------------------------------- histograms --- */

/* Count occurrences of each byte value 0-255. */
void count_bytes(const uint8_t *data, size_t len, long freq[256])
{
    memset(freq, 0, 256 * sizeof freq[0]);
    for (size_t i = 0; i < len; i++) {
        freq[data[i]]++;
    }
}

/*
 * Shannon entropy H = -Σ p(x)·log2(p(x)) in bits per byte, from a frequency
 * histogram. Zero-count bins contribute nothing. `total` must be the sum of
 * `frequencies`.
 */
double shannon_bits_per_byte(const long *frequencies, size_t bins, double total)
{
    double h = 0.0;

    if (total <= 0.0) {
        return 0.0;
    }
    for (size_t i = 0; i < bins; i++) {
        double p;

        if (frequencies[i] == 0) {
            continue;
        }
        p = (double) frequencies[i] / total;
        h -= p * log2(p);
    }
    return h;
}

/*
 * Pearson χ² comparing observed byte counts against a uniform expectation
 * E = total/256 per bin. Returns NAN for a non-positive sample size (the TS
 * reference throws; a C caller checks with isnan()).
 */
double chi_squared_statistic(const long *frequencies, double total)
{
    double expected, chi2 = 0.0;

    if (total <= 0.0) {
        return NAN; /* chi-squared needs a positive sample size */
    }
    expected = total / 256.0;
    for (size_t i = 0; i < 256; i++) {
        double diff = (double) frequencies[i] - expected;
        chi2 += (diff * diff) / expected;
    }
    return chi2;
}

/* --------------------------------------------------------- gamma functions --- */

/* Lanczos approximation (g=7, 9 coefficients) of ln Γ(x). */
static double ln_gamma(double x)
{
    static const double g[9] = {
        0.99999999999980993,  676.5203681218851,   -1259.1392167224028,
        771.32342877765313,  -176.61502916214059,     12.507343278686905,
        -0.13857109526572012,   9.9843695780195716e-6, 1.5056327351493116e-7
    };
    double a, t;

    if (x < 0.5) {
        /* Reflection formula: Γ(x)·Γ(1-x) = π / sin(πx) */
        return log(M_PI / sin(M_PI * x)) - ln_gamma(1.0 - x);
    }
    x -= 1.0;
    a = g[0];
    t = x + 7.5;
    for (int i = 1; i < 9; i++) {
        a += g[i] / (x + i);
    }
    return 0.5 * log(2.0 * M_PI) + (x + 0.5) * log(t) - t + log(a);
}

static double clamp01(double v)
{
    return (v < 0.0) ? 0.0 : (v > 1.0) ? 1.0 : v;
}

/*
 * Regularized upper incomplete gamma function Q(a, x) = Γ(a,x)/Γ(a), via the
 * power series (x < a+1) or the Lentz continued fraction (otherwise).
 * Numerical Recipes §6.2.
 */
double gamma_q(double a, double x)
{
    if (a <= 0.0 || x < 0.0) {
        return NAN;
    }
    if (x == 0.0) {
        return 1.0;
    }

    if (x < a + 1.0) {
        /* Series for P(a,x); Q = 1 - P */
        double ap = a;
        double sum = 1.0 / a;
        double del = sum;

        for (int n = 0; n < 1000; n++) {
            ap += 1.0;
            del *= x / ap;
            sum += del;
            if (fabs(del) < fabs(sum) * 1e-15) {
                break;
            }
        }
        return clamp01(1.0 - sum * exp(-x + a * log(x) - ln_gamma(a)));
    }

    /* Continued fraction for Q(a,x) */
    {
        const double FPMIN = 1e-300;
        double b = x + 1.0 - a;
        double c = 1.0 / FPMIN;
        double d = 1.0 / b;
        double h = d;

        for (int i = 1; i <= 1000; i++) {
            double an = -i * (i - a);
            double del;

            b += 2.0;
            d = an * d + b;
            if (fabs(d) < FPMIN) d = FPMIN;
            c = b + an / c;
            if (fabs(c) < FPMIN) c = FPMIN;
            d = 1.0 / d;
            del = d * c;
            h *= del;
            if (fabs(del - 1.0) < 1e-15) {
                break;
            }
        }
        return clamp01(exp(-x + a * log(x) - ln_gamma(a)) * h);
    }
}

/* Upper-tail p-value for a χ² statistic with `df` degrees of freedom. */
double chi_squared_p(double chi2, double df)
{
    if (chi2 < 0.0 || df <= 0.0) {
        return NAN;
    }
    return gamma_q(df / 2.0, chi2 / 2.0);
}

/* ---------------------------------------------------------------- measures --- */

/*
 * Circular serial correlation between consecutive bytes (the `ent` tool's
 * metric): scc = (Σxy - (Σx)²/n) / (Σx² - (Σx)²/n) over the pair sequence
 * (x₀,x₁), (x₁,x₂), …, (xₙ₋₁,x₀). 0 = uncorrelated, ±1 = perfectly
 * (anti)correlated. Constant input has a zero denominator → reported as 0
 * (nothing to correlate); inputs shorter than 2 bytes are also 0.
 */
double serial_correlation_coefficient(const uint8_t *data, size_t n)
{
    double sum = 0.0, sum_sq = 0.0, sum_xy = 0.0, mean_sq, denom;

    if (n < 2) {
        return 0.0;
    }
    for (size_t i = 0; i < n; i++) {
        double x = data[i];
        double y = data[(i + 1) % n];

        sum    += x;
        sum_sq += x * x;
        sum_xy += x * y;
    }
    mean_sq = (sum * sum) / (double) n;
    denom = sum_sq - mean_sq;
    if (denom == 0.0) {
        return 0.0;
    }
    return (sum_xy - mean_sq) / denom;
}

/*
 * Monte Carlo π estimate: consecutive byte pairs are (x, y) points in a
 * 256×256 square; the fraction inside the inscribed circle (center 127.5,
 * radius 128) times 4 estimates π. Inputs with fewer than 2 bytes → 0.
 */
double monte_carlo_pi_estimate(const uint8_t *data, size_t len)
{
    size_t pairs = len / 2;
    long inside = 0;

    if (pairs == 0) {
        return 0.0;
    }
    for (size_t i = 0; i < pairs; i++) {
        double dx = (double) data[2 * i]     - 127.5;
        double dy = (double) data[2 * i + 1] - 127.5;

        if (dx * dx + dy * dy <= 128.0 * 128.0) {
            inside++;
        }
    }
    return (4.0 * (double) inside) / (double) pairs;
}

/*
 * Shannon entropy (bits/byte) of each consecutive `block_size`-byte block.
 * Returns a malloc'd array of length *out_count, or NULL when block_size < 1.
 */
double *block_entropies(const uint8_t *data, size_t len, size_t block_size, size_t *out_count)
{
    double *blocks;
    size_t count;

    if (block_size < 1) {
        *out_count = 0;
        return NULL; /* blockSize must be at least 1 */
    }
    count = (len + block_size - 1) / block_size;
    *out_count = count;
    if (count == 0) {
        return NULL;
    }
    blocks = malloc(count * sizeof *blocks);
    if (blocks == NULL) {
        *out_count = 0;
        return NULL;
    }

    for (size_t b = 0; b < count; b++) {
        long counts[256];
        size_t off = b * block_size;
        size_t end = off + block_size;
        size_t n;

        if (end > len) {
            end = len;
        }
        n = end - off;
        count_bytes(data + off, n, counts);
        blocks[b] = shannon_bits_per_byte(counts, 256, (double) n);
    }
    return blocks;
}

/* Map a Shannon entropy (bits/byte) to the verdict scale. */
entropy_verdict verdict_from_shannon(double bits_per_byte)
{
    if (bits_per_byte > 7.5) return VERDICT_EXCELLENT;
    if (bits_per_byte > 6.0) return VERDICT_GOOD;
    if (bits_per_byte > 4.0) return VERDICT_SUSPICIOUS;
    return VERDICT_LOW;
}

/* ---------------------------------------------------------------- analysis --- */

/*
 * Full analysis of a byte sequence. Returns false on empty input — there is
 * nothing to measure and every metric would be undefined. On success the caller
 * owns `out->block_entropies` and must free_entropy_analysis() it.
 */
bool analyze_entropy(const uint8_t *data, size_t len, entropy_analysis *out)
{
    if (data == NULL || len == 0 || out == NULL) {
        return false; /* Nothing to analyze - provide at least 1 byte of data. */
    }
    memset(out, 0, sizeof *out);

    count_bytes(data, len, out->byte_frequencies);
    out->shannon_entropy     = shannon_bits_per_byte(out->byte_frequencies, 256, (double) len);
    out->chi_squared         = chi_squared_statistic(out->byte_frequencies, (double) len);
    out->chi_squared_p_value = chi_squared_p(out->chi_squared, CHI_SQUARED_DF);
    out->serial_correlation  = serial_correlation_coefficient(data, len);
    out->monte_carlo_pi      = monte_carlo_pi_estimate(data, len);
    out->block_entropies     = block_entropies(data, len, BLOCK_SIZE, &out->block_count);
    out->verdict             = verdict_from_shannon(out->shannon_entropy);
    return true;
}

void free_entropy_analysis(entropy_analysis *analysis)
{
    if (analysis == NULL) {
        return;
    }
    free(analysis->block_entropies);
    analysis->block_entropies = NULL;
    analysis->block_count = 0;
}

/* -------------------------------------------------------------------- demo --- */

static void report(const char *label, const uint8_t *data, size_t len)
{
    entropy_analysis a;

    if (!analyze_entropy(data, len, &a)) {
        printf("%-18s (empty - nothing to analyze)\n", label);
        return;
    }
    printf("%-18s shannon %.4f bits/byte  chi2 %9.1f  p %.4f  scc %+.4f  pi %.4f  -> %s\n",
           label, a.shannon_entropy, a.chi_squared, a.chi_squared_p_value,
           a.serial_correlation, a.monte_carlo_pi, verdict_name(a.verdict));
    free_entropy_analysis(&a);
}

int main(void)
{
    uint8_t uniform[4096];
    uint8_t constant[4096];
    uint8_t ascii[4096];

    /* Every byte value equally often: maximal entropy, χ² of exactly 0. */
    for (size_t i = 0; i < sizeof uniform; i++) {
        uniform[i] = (uint8_t) (i % 256);
    }
    /* One repeating byte: zero entropy. */
    memset(constant, 0x41, sizeof constant);
    /* Printable ASCII only: ~6.6 bits/byte — "good" by Shannon alone, yet the
     * χ² p-value exposes it immediately. This is the case the multi-metric
     * battery exists to catch. */
    for (size_t i = 0; i < sizeof ascii; i++) {
        ascii[i] = (uint8_t) (32 + (i * 7 + i / 95) % 95);
    }

    report("uniform ramp:",  uniform,  sizeof uniform);
    report("constant 0x41:", constant, sizeof constant);
    report("printable ASCII:", ascii,  sizeof ascii);

    /* The verdict scale, for reference. */
    printf("\nverdict thresholds (bits/byte): >7.5 excellent, >6.0 good, >4.0 suspicious, else low\n");
    return EXIT_SUCCESS;
}