Entropy Visualizer — Zig source

//! entropy-visualizer — byte-level randomness analysis. No dependencies.
//!
//! Language: Zig 0.14 (standard library only)
//! Ported from: src/lib/entropy-visualizer.ts (the canonical TypeScript implementation).
//! display source — part of CosmoDev's polyglot tool pages.
//!
//! Everything is deterministic: the same bytes always produce the same
//! numbers. (This file is the contract the Go CLI twin mirrors — see CLAUDE.md
//! "Dual source".)

const std = @import("std");

pub const EntropyVerdict = enum {
    excellent,
    good,
    suspicious,
    low,
};

/// Bytes per heatmap block.
pub const BLOCK_SIZE = 16;

/// Degrees of freedom for the byte-frequency χ² test (256 bins - 1 constraint).
pub const CHI_SQUARED_DF = 255;

pub const EntropyAnalysis = struct {
    /// Shannon entropy in bits per byte (0 = one repeating byte, 8 = perfectly uniform).
    shannon_entropy: f64,
    /// χ² statistic against the uniform 256-bin expectation.
    chi_squared: f64,
    /// Upper-tail p-value for χ² with 255 degrees of freedom (1 = perfectly plausible).
    chi_squared_p_value: f64,
    /// Circular serial correlation between consecutive bytes (-1 … +1, 0 = uncorrelated).
    serial_correlation: f64,
    /// Monte Carlo π estimate from consecutive byte pairs (≈3.14159 for random data).
    monte_carlo_pi: f64,
    /// Occurrence count per byte value 0-255 (always 256 entries).
    byte_frequencies: [256]u64,
    /// Shannon entropy (bits/byte) of each BLOCK_SIZE-byte block, for the heatmap.
    block_entropies: []f64,
    verdict: EntropyVerdict,
};

pub const Error = error{
    EmptyInput,
    InvalidBlockSize,
    NonPositiveSample,
    OutOfMemory,
};

/// Count occurrences of each byte value 0-255.
pub fn countBytes(freq: *[256]u64, data: []const u8) void {
    @memset(freq, 0);
    for (data) |b| freq[b] += 1;
}

/// Shannon entropy H = -Σ p(x)·log2(p(x)) in bits per byte, computed from a
/// frequency histogram. Zero-count bins contribute nothing. `total` must be the
/// sum of `frequencies`.
pub fn shannonBitsPerByte(frequencies: *const [256]u64, total: u64) f64 {
    if (total == 0) return 0;
    const t = @as(f64, @floatFromInt(total));
    var h: f64 = 0;
    for (frequencies) |count| {
        if (count == 0) continue;
        const p = @as(f64, @floatFromInt(count)) / t;
        h -= p * std.math.log2(p);
    }
    return h;
}

/// Pearson χ² comparing observed byte counts against a uniform expectation
/// E = total/256 per bin.
pub fn chiSquaredStatistic(frequencies: *const [256]u64, total: u64) Error!f64 {
    if (total == 0) return Error.NonPositiveSample;
    const expected = @as(f64, @floatFromInt(total)) / 256.0;
    var chi2: f64 = 0;
    for (frequencies) |observed| {
        const diff = @as(f64, @floatFromInt(observed)) - expected;
        chi2 += (diff * diff) / expected;
    }
    return chi2;
}

/// Lanczos approximation (g=7, 9 coefficients) of ln Γ(x).
fn lnGamma(x_in: f64) f64 {
    const g = [9]f64{
        0.99999999999980993,  676.5203681218851,   -1259.1392167224028,
        771.32342877765313,   -176.61502916214059, 12.507343278686905,
        -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
    };
    var x = x_in;
    if (x < 0.5) {
        // Reflection formula: Γ(x)·Γ(1-x) = π / sin(πx)
        return @log(std.math.pi / @sin(std.math.pi * x)) - lnGamma(1 - x);
    }
    x -= 1;
    var a = g[0];
    const t = x + 7.5;
    for (1..9) |i| {
        a += g[i] / (x + @as(f64, @floatFromInt(i)));
    }
    return 0.5 * @log(2 * std.math.pi) + (x + 0.5) * @log(t) - t + @log(a);
}

/// 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.
pub fn gammaQ(a: f64, x: f64) f64 {
    if (a <= 0 or x < 0) return std.math.nan(f64);
    if (x == 0) return 1;
    if (x < a + 1) {
        // Series for P(a,x); Q = 1 - P
        var ap = a;
        var sum = 1 / a;
        var del = sum;
        for (0..1000) |_| {
            ap += 1;
            del *= x / ap;
            sum += del;
            if (@abs(del) < @abs(sum) * 1e-15) break;
        }
        return @min(1, @max(0, 1 - sum * @exp(-x + a * @log(x) - lnGamma(a))));
    }
    // Continued fraction for Q(a,x)
    const FPMIN = 1e-300;
    var b = x + 1 - a;
    var c = 1 / FPMIN;
    var d = 1 / b;
    var h = d;
    for (1..1001) |i| {
        const fi: f64 = @floatFromInt(i);
        const an = -fi * (fi - a);
        b += 2;
        d = an * d + b;
        if (@abs(d) < FPMIN) d = FPMIN;
        c = b + an / c;
        if (@abs(c) < FPMIN) c = FPMIN;
        d = 1 / d;
        const del = d * c;
        h *= del;
        if (@abs(del - 1) < 1e-15) break;
    }
    return @min(1, @max(0, @exp(-x + a * @log(x) - lnGamma(a)) * h));
}

/// Upper-tail p-value for a χ² statistic with `df` degrees of freedom.
pub fn chiSquaredP(chi2: f64, df: f64) f64 {
    if (chi2 < 0 or df <= 0) return std.math.nan(f64);
    return gammaQ(df / 2, chi2 / 2);
}

/// 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.
pub fn serialCorrelationCoefficient(data: []const u8) f64 {
    const n = data.len;
    if (n < 2) return 0;
    const fn_n: f64 = @floatFromInt(n);
    var sum: f64 = 0;
    var sum_sq: f64 = 0;
    var sum_xy: f64 = 0;
    for (data, 0..) |x, i| {
        const y = data[(i + 1) % n];
        const fx: f64 = @floatFromInt(x);
        const fy: f64 = @floatFromInt(y);
        sum += fx;
        sum_sq += fx * fx;
        sum_xy += fx * fy;
    }
    const mean_sq = (sum * sum) / fn_n;
    const denom = sum_sq - mean_sq;
    if (denom == 0) return 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.
pub fn monteCarloPiEstimate(data: []const u8) f64 {
    const pairs = data.len / 2;
    if (pairs == 0) return 0;
    var inside: f64 = 0;
    for (0..pairs) |i| {
        const dx = @as(f64, @floatFromInt(data[2 * i])) - 127.5;
        const dy = @as(f64, @floatFromInt(data[2 * i + 1])) - 127.5;
        if (dx * dx + dy * dy <= 128 * 128) inside += 1;
    }
    return (4 * inside) / @as(f64, @floatFromInt(pairs));
}

/// Shannon entropy (bits/byte) of each consecutive `block_size`-byte block.
/// Caller owns the returned slice.
pub fn blockEntropies(allocator: std.mem.Allocator, data: []const u8, block_size: usize) Error![]f64 {
    if (block_size < 1) return Error.InvalidBlockSize;
    var blocks = std.ArrayList(f64).init(allocator);
    errdefer blocks.deinit();
    var counts: [256]u64 = undefined;
    var off: usize = 0;
    while (off < data.len) : (off += block_size) {
        @memset(&counts, 0);
        const end = @min(off + block_size, data.len);
        for (data[off..end]) |b| counts[b] += 1;
        try blocks.append(shannonBitsPerByte(&counts, end - off));
    }
    return blocks.toOwnedSlice();
}

/// Map a Shannon entropy (bits/byte) to the verdict scale.
pub fn verdictFromShannon(bits_per_byte: f64) EntropyVerdict {
    if (bits_per_byte > 7.5) return .excellent;
    if (bits_per_byte > 6.0) return .good;
    if (bits_per_byte > 4.0) return .suspicious;
    return .low;
}

/// Full analysis of a byte sequence. Errors on empty input — there is nothing
/// to measure and every metric would be undefined. Caller owns
/// `block_entropies` (free with the same allocator).
pub fn analyzeEntropy(allocator: std.mem.Allocator, data: []const u8) Error!EntropyAnalysis {
    if (data.len == 0) return Error.EmptyInput;
    const total: u64 = data.len;
    var freq: [256]u64 = undefined;
    countBytes(&freq, data);
    const shannon = shannonBitsPerByte(&freq, total);
    const chi2 = try chiSquaredStatistic(&freq, total);
    return .{
        .shannon_entropy = shannon,
        .chi_squared = chi2,
        .chi_squared_p_value = chiSquaredP(chi2, CHI_SQUARED_DF),
        .serial_correlation = serialCorrelationCoefficient(data),
        .monte_carlo_pi = monteCarloPiEstimate(data),
        .byte_frequencies = freq,
        .block_entropies = try blockEntropies(allocator, data, BLOCK_SIZE),
        .verdict = verdictFromShannon(shannon),
    };
}