src/n_entropy.c

Implementation of the byte-sample randomness metrics.

Source code

/*
 * Nilorea Library
 * Copyright (C) 2005-2026 Castagnier Mickael
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
 * implied. See the License for the specific language governing
 * permissions and limitations under the License.
 *
 * SPDX-License-Identifier: Apache-2.0
 */


#include "nilorea/n_entropy.h"

#include <math.h>

double n_entropy_shannon(const unsigned char* data, size_t len) {
    size_t counts[256];
    size_t i;
    double h = 0.0;
    if (!data || len == 0)
        return 0.0;
    for (i = 0; i < 256; i++)
        counts[i] = 0;
    for (i = 0; i < len; i++)
        counts[data[i]]++;
    for (i = 0; i < 256; i++) {
        if (counts[i]) {
            double p = (double)counts[i] / (double)len;
            h -= p * log2(p);
        }
    }
    return h;
}

double n_entropy_monobit(const unsigned char* data, size_t len) {
    size_t i;
    size_t ones = 0;
    if (!data || len == 0)
        return 0.0;
    for (i = 0; i < len; i++) {
        unsigned int b = data[i];
        /* popcount of one byte */
        b = b - ((b >> 1) & 0x55u);
        b = (b & 0x33u) + ((b >> 2) & 0x33u);
        b = (b + (b >> 4)) & 0x0Fu;
        ones += b;
    }
    return (double)ones / ((double)len * 8.0);
}

double n_entropy_chi_square(const unsigned char* data, size_t len) {
    size_t counts[256];
    size_t i;
    double expected, chi = 0.0;
    if (!data || len == 0)
        return 0.0;
    for (i = 0; i < 256; i++)
        counts[i] = 0;
    for (i = 0; i < len; i++)
        counts[data[i]]++;
    expected = (double)len / 256.0;
    for (i = 0; i < 256; i++) {
        double diff = (double)counts[i] - expected;
        chi += (diff * diff) / expected;
    }
    return chi;
}

double n_entropy_chi_square_p(double chi, size_t df) {
    double d, t, z;
    if (df == 0 || chi <= 0.0)
        return 1.0;
    d = (double)df;
    /* Wilson-Hilferty: (chi2/df)^(1/3) is close to normal with mean 1-2/(9df)
       and variance 2/(9df), so the upper tail becomes a normal tail */
    t = cbrt(chi / d);
    z = (t - (1.0 - 2.0 / (9.0 * d))) / sqrt(2.0 / (9.0 * d));
    return 0.5 * erfc(z / sqrt(2.0));
}

size_t n_entropy_positions(const unsigned char* const* samples, size_t count, size_t sample_len, N_ENTROPY_POS* out, size_t max_out) {
    size_t pos, n;
    if (!samples || !out || count == 0 || sample_len == 0 || max_out == 0)
        return 0;
    n = (sample_len < max_out) ? sample_len : max_out;
    for (pos = 0; pos < n; pos++) {
        size_t counts[256];
        size_t i, distinct = 0;
        double h = 0.0, chi = 0.0;
        for (i = 0; i < 256; i++)
            counts[i] = 0;
        for (i = 0; i < count; i++)
            if (samples[i])
                counts[samples[i][pos]]++;
        for (i = 0; i < 256; i++)
            if (counts[i]) {
                double p = (double)counts[i] / (double)count;
                distinct++;
                h -= p * log2(p);
            }
        /* the chi-square is taken over the alphabet ACTUALLY seen, not over all 256
           values: a hex token uses 16 of them by design, and scoring it against 256
           would report every well-formed token as catastrophically biased */
        if (distinct > 0) {
            double expected = (double)count / (double)distinct;
            for (i = 0; i < 256; i++)
                if (counts[i]) {
                    double diff = (double)counts[i] - expected;
                    chi += (diff * diff) / expected;
                }
        }
        out[pos].distinct = distinct;
        out[pos].shannon = h;
        out[pos].chi_square = chi;
        out[pos].p_value = n_entropy_chi_square_p(chi, distinct > 0 ? distinct - 1 : 0);
    }
    return n;
}

size_t n_entropy_bits(const unsigned char* const* samples, size_t count, size_t sample_len, N_ENTROPY_BIT* out, size_t max_out) {
    size_t bit, nbits;
    if (!samples || !out || count == 0 || sample_len == 0 || max_out == 0)
        return 0;
    nbits = sample_len * 8;
    if (nbits > max_out)
        nbits = max_out;
    for (bit = 0; bit < nbits; bit++) {
        size_t byte = bit / 8;
        unsigned mask = (unsigned)(0x80u >> (bit % 8));
        size_t i, ones = 0;
        double frac, z;
        for (i = 0; i < count; i++)
            if (samples[i] && (samples[i][byte] & mask))
                ones++;
        frac = (double)ones / (double)count;
        /* z of the imbalance against a fair coin: (ones - n/2) / (sqrt(n)/2) */
        z = ((double)ones - (double)count / 2.0) / (sqrt((double)count) / 2.0);
        out[bit].ones = ones;
        out[bit].fraction = frac;
        out[bit].z = z;
        out[bit].p_value = erfc(fabs(z) / sqrt(2.0));
    }
    return nbits;
}

double n_entropy_runs_p(const unsigned char* data, size_t len) {
    size_t nbits, i, runs = 1;
    double pi, z;
    int prev;
    if (!data || len == 0)
        return 0.0;
    nbits = len * 8;
    pi = n_entropy_monobit(data, len);
    /* the test is only meaningful when the stream is roughly balanced; outside that
       band monobit has already failed and this would divide by a vanishing term */
    if (fabs(pi - 0.5) >= 2.0 / sqrt((double)nbits))
        return 0.0;
    prev = (data[0] & 0x80u) ? 1 : 0;
    for (i = 1; i < nbits; i++) {
        int b = (data[i / 8] & (unsigned)(0x80u >> (i % 8))) ? 1 : 0;
        if (b != prev)
            runs++;
        prev = b;
    }
    z = fabs((double)runs - 2.0 * (double)nbits * pi * (1.0 - pi)) /
        (2.0 * sqrt(2.0 * (double)nbits) * pi * (1.0 - pi));
    return erfc(z / sqrt(2.0));
}

size_t n_entropy_longest_run(const unsigned char* data, size_t len) {
    size_t nbits, i, best = 0, cur = 1;
    int prev;
    if (!data || len == 0)
        return 0;
    nbits = len * 8;
    prev = (data[0] & 0x80u) ? 1 : 0;
    best = 1;
    for (i = 1; i < nbits; i++) {
        int b = (data[i / 8] & (unsigned)(0x80u >> (i % 8))) ? 1 : 0;
        if (b == prev) {
            cur++;
            if (cur > best)
                best = cur;
        } else {
            cur = 1;
        }
        prev = b;
    }
    return best;
}