#!/home/ben/software/install/bin/perl
use warnings;
use strict;

# A small error factor.

use constant {
    eps => 1e-10,
};

print "Entropy of a single constant:\n";

print &entropy (1), "\n\n";

print "Entropy of alphabet with probabilities 1/3/, 1/4, 1/4, 1/6:\n";

my @p = (1/3, 1/4, 1/4, 1/6);

print &entropy (@p), "\n\n";

print "Entropy of four equally-likely symbols:\n";

my @p3 = (1/4) x 4;

print &entropy (@p3), "\n\n";

print "Entropy of six equally-likely symbols:\n";

my @p2 = (1/6) x 6;

print &entropy (@p2), " (log_2 (6) = ", log2 (6), ")\n\n";

print "Entropy of eight equally-likely symbols:\n";

my @p4 = (1/8) x 8;

print &entropy (@p4), "\n\n";

print "Entropy of a two-symbol alphabet with probabilities 1/n, 1 - 1/n:\n";

for my $factor (1..10) {
    my @p5 = (1 - 1/$factor, 1/$factor);
    print $factor, " ", entropy (@p5), "\n";
}

print "\n";

print "Entropy of an n-symbol alphabet, probs (1/2, 1/4, ... 1/(2^n), 1/(2^n)):\n";

for my $factor (1..10) {
    my @q;
    for my $j (0..$factor) {
	$q[$j] = 1/(2**($j+1));
    }
    $q[$factor + 1] = 1/(2**($factor+1));
    printf "n = %d: entropy = %g\n", $factor + 1, entropy (@q);
}

print "\n";

exit;

# Log to the base 2, so log2 (2) = 1, log2 (4) = 2, etc.

sub log2
{
    my ($n) = @_;
    return log ($n) / log (2);
}

# Given a set of discrete probabilities, "@p", compute the entropy.

sub entropy
{
    my (@p) = @_;

    # Do some checks on the inputs before proceeding.

    my $sum = 0;
    for (@p) {
	if ($_ < 0 || $_ > 1) {
	    warn "$_ is not a probability.\n";
	    return undef;
	}
	$sum += $_;
    }
    if (abs ($sum - 1) > eps) {
	warn "@p does not sum to 1.\n";
	return undef;
    }

    # The return value.

    my $h = 0;

    for (@p) {
	if ($_ == 0) {
	    next;
	}
	$h += - $_ * log2 ($_);
    }
    return $h;
}