Entropy Visualizer — C++ source

// Entropy Visualizer — byte-level randomness analysis.
//
// Language: C++17 (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.
// Shannon entropy, Pearson χ² against the uniform expectation (with the
// upper-tail p-value via the regularized incomplete gamma function, Numerical
// Recipes §6.2), the `ent` tool's circular serial correlation, and a Monte
// Carlo π estimate — the same battery `ent` runs on a file.

#include <array>
#include <cmath>
#include <cstdint>
#include <stdexcept>
#include <string>
#include <vector>

namespace entropy {

using Bytes = std::vector<uint8_t>;
using Frequencies = std::array<int, 256>;

enum class Verdict { Excellent, Good, Suspicious, Low };

/** Bytes per heatmap block. */
constexpr int BLOCK_SIZE = 16;
/** Degrees of freedom for the byte-frequency χ² test (256 bins − 1 constraint). */
constexpr int CHI_SQUARED_DF = 255;

struct EntropyAnalysis {
  double shannonEntropy;       ///< bits per byte (0 = one repeating byte, 8 = uniform)
  double chiSquared;           ///< χ² statistic against the uniform 256-bin expectation
  double chiSquaredPValue;     ///< upper-tail p-value, χ² with 255 df (1 = perfectly plausible)
  double serialCorrelation;    ///< circular correlation between consecutive bytes (−1 … +1)
  double monteCarloPi;         ///< π estimate from consecutive byte pairs (≈3.14159 for random data)
  Frequencies byteFrequencies; ///< occurrence count per byte value 0-255
  std::vector<double> blockEntropies; ///< per-block Shannon entropy, for the heatmap
  Verdict verdict;
};

/** Count occurrences of each byte value 0-255. */
Frequencies countBytes(const Bytes& data) {
  Frequencies freq{};
  for (uint8_t b : data) freq[b]++;
  return freq;
}

/**
 * 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 shannonBitsPerByte(const Frequencies& frequencies, long long total) {
  if (total <= 0) return 0;
  double h = 0;
  for (int count : frequencies) {
    if (count == 0) continue;
    const double p = static_cast<double>(count) / static_cast<double>(total);
    h -= p * std::log2(p);
  }
  return h;
}

/** Pearson χ² comparing observed byte counts against a uniform expectation
 *  E = total/256 per bin. */
double chiSquaredStatistic(const Frequencies& frequencies, long long total) {
  if (total <= 0) throw std::invalid_argument("chi-squared needs a positive sample size");
  const double expected = static_cast<double>(total) / 256.0;
  double chi2 = 0;
  for (int observed : frequencies) {
    const double diff = static_cast<double>(observed) - expected;
    chi2 += (diff * diff) / expected;
  }
  return chi2;
}

/** Lanczos approximation (g=7, 9 coefficients) of ln Γ(x). */
double lnGamma(double x) {
  static const double g[] = {
    0.99999999999980993, 676.5203681218851, -1259.1392167224028,
    771.32342877765313, -176.61502916214059, 12.507343278686905,
    -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
  };
  if (x < 0.5) {
    // Reflection formula: Γ(x)·Γ(1−x) = π / sin(πx)
    return std::log(M_PI / std::sin(M_PI * x)) - lnGamma(1 - x);
  }
  x -= 1;
  double a = g[0];
  const double t = x + 7.5;
  for (int i = 1; i < 9; i++) a += g[i] / (x + i);
  return 0.5 * std::log(2 * M_PI) + (x + 0.5) * std::log(t) - t + std::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.
 */
double gammaQ(double a, double x) {
  const double NAN_ = std::nan("");
  if (a <= 0 || x < 0) return NAN_;
  if (x == 0) return 1;
  if (x < a + 1) {
    // Series for P(a,x); Q = 1 − P
    double ap = a;
    double sum = 1 / a;
    double del = sum;
    for (int n = 0; n < 1000; n++) {
      ap += 1;
      del *= x / ap;
      sum += del;
      if (std::fabs(del) < std::fabs(sum) * 1e-15) break;
    }
    return std::min(1.0, std::max(0.0, 1 - sum * std::exp(-x + a * std::log(x) - lnGamma(a))));
  }
  // Continued fraction for Q(a,x)
  constexpr double FPMIN = 1e-300;
  double b = x + 1 - a;
  double c = 1 / FPMIN;
  double d = 1 / b;
  double h = d;
  for (int i = 1; i <= 1000; i++) {
    const double an = -static_cast<double>(i) * (i - a);
    b += 2;
    d = an * d + b;
    if (std::fabs(d) < FPMIN) d = FPMIN;
    c = b + an / c;
    if (std::fabs(c) < FPMIN) c = FPMIN;
    d = 1 / d;
    const double del = d * c;
    h *= del;
    if (std::fabs(del - 1) < 1e-15) break;
  }
  return std::min(1.0, std::max(0.0, std::exp(-x + a * std::log(x) - lnGamma(a)) * h));
}

/** Upper-tail p-value for a χ² statistic with `df` degrees of freedom. */
double chiSquaredP(double chi2, int df = CHI_SQUARED_DF) {
  if (chi2 < 0 || df <= 0) return std::nan("");
  return gammaQ(df / 2.0, chi2 / 2.0);
}

/**
 * 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;
 * inputs shorter than 2 bytes are also 0.
 */
double serialCorrelationCoefficient(const Bytes& data) {
  const size_t n = data.size();
  if (n < 2) return 0;
  long long sum = 0, sumSq = 0, sumXY = 0;
  for (size_t i = 0; i < n; i++) {
    const long long x = data[i];
    const long long y = data[(i + 1) % n];
    sum += x;
    sumSq += x * x;
    sumXY += x * y;
  }
  const double meanSq = static_cast<double>(sum * sum) / static_cast<double>(n);
  const double denom = static_cast<double>(sumSq) - meanSq;
  if (denom == 0) return 0;
  return (static_cast<double>(sumXY) - meanSq) / 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 monteCarloPiEstimate(const Bytes& data) {
  const size_t pairs = data.size() / 2;
  if (pairs == 0) return 0;
  long long inside = 0;
  for (size_t i = 0; i < pairs; i++) {
    const double dx = static_cast<double>(data[2 * i]) - 127.5;
    const double dy = static_cast<double>(data[2 * i + 1]) - 127.5;
    if (dx * dx + dy * dy <= 128.0 * 128.0) inside++;
  }
  return (4.0 * static_cast<double>(inside)) / static_cast<double>(pairs);
}

/** Shannon entropy (bits/byte) of each consecutive `blockSize`-byte block. */
std::vector<double> blockEntropies(const Bytes& data, int blockSize = BLOCK_SIZE) {
  if (blockSize < 1) throw std::invalid_argument("blockSize must be at least 1");
  std::vector<double> blocks;
  for (size_t off = 0; off < data.size(); off += static_cast<size_t>(blockSize)) {
    Frequencies counts{};
    long long n = 0;
    const size_t end = std::min(off + static_cast<size_t>(blockSize), data.size());
    for (size_t i = off; i < end; i++) {
      counts[data[i]]++;
      n++;
    }
    blocks.push_back(shannonBitsPerByte(counts, n));
  }
  return blocks;
}

/** Map a Shannon entropy (bits/byte) to the verdict scale. */
Verdict verdictFromShannon(double bitsPerByte) {
  if (bitsPerByte > 7.5) return Verdict::Excellent;
  if (bitsPerByte > 6.0) return Verdict::Good;
  if (bitsPerByte > 4.0) return Verdict::Suspicious;
  return Verdict::Low;
}

/**
 * Full analysis of a byte sequence. Throws on empty input — there is nothing
 * to measure and every metric would be undefined.
 */
EntropyAnalysis analyzeEntropy(const Bytes& data) {
  if (data.empty()) {
    throw std::invalid_argument("Nothing to analyze - provide at least 1 byte of data.");
  }
  const long long total = static_cast<long long>(data.size());
  const Frequencies freq = countBytes(data);
  const double shannon = shannonBitsPerByte(freq, total);
  const double chi2 = chiSquaredStatistic(freq, total);
  return EntropyAnalysis{
    shannon,
    chi2,
    chiSquaredP(chi2),
    serialCorrelationCoefficient(data),
    monteCarloPiEstimate(data),
    freq,
    blockEntropies(data),
    verdictFromShannon(shannon),
  };
}

} // namespace entropy