#!/usr/local/bin/perl -w
#hz.p program to find the entropy of z when it is given by A x = z 
# using GF(q) for q in 2,4,8,16,32,64
#
# To avoid an exhaustive search, we use random sampling.

my $usage = "$0 [-q q] [-r row-weight] [-f bsc-noise-level] [-N iterations] [-v]\n\
\tdefaults: q=8, r=6, f=0.008, N=5000, verbosity off.\n";

$tr = 6 ; # number of elemnts t  oadd.
$f = 0.08 ; # noise of BSC
$q =   8 ; 
$init1 = 0; # Do we set first element to 1, always?
$verbose = 0 ;
$exhaustive = 1 ; 
$qary = 0 ;
$Runs=5000;

my $option;
while (defined($ARGV[0]) && $ARGV[0] =~ /^-/) {
	$option = shift();
	if ($option eq "-q") {
		$q = shift();
		die "unsupported field size: $q.  Use 2, 4, 8, ... 64.\n"
			unless (($q == 2) || ($q == 4) || ($q == 8) 
				|| ($q == 16) || ($q == 32) || ($q == 64));
	} elsif ($option eq "-r") {
		$tr = shift();
		die "row weight must be at least 2\n" unless
			($tr >= 2);
		
	} elsif ($option eq "-f") {
		$f = shift();
		die "BSC noise level must be in range 0 < f < 1\n" unless
			(($f > 0) && ($f < 1));
	} elsif ($option eq "-N") {
		$Runs = shift();
		die $usage unless ($Runs > 0);
	} elsif ($option eq "-b") {
		$verbose = 1;
	} else {
		die $usage;
	}
}

$Q = $q-1 ; 
@LOG4 = (0,0,1,2);
@LOG8 = (0,0,1,3,2,6,4,5);
@LOG16= (0,0,1,4,2,8,5,10,3,14,9,7,6,13,11,12);
@LOG32= (0,0,1,18,2,5,19,11,3,29,6,27,20,8,12,23,4,
	 10,30,17,7,22,28,26,21,25,9,16,13,14,24,15);
@LOG64= (0,0,1,6,2,12,7,26,3,32,13,35,8,48,27,18,4,24,
	 33,16,14,52,36,54,9,45,49,38,28,41,19,56,5,62,
	 25,11,34,31,17,47,15,23,53,51,37,44,55,40,10,
	 61,46,30,50,22,39,43,29,60,42,21,20,59,57,58);

@EXP4 = (1,2,3);
@EXP8 = (1,2,4,3,6,7,5);
@EXP16= (1,2,4,8,3,6,12,11,5,10,7,14,15,13,9);
@EXP32= (1,2,4,8,16,5,10,20,13,26,17,7,14,28,29,31,
	 27,19,3,6,12,24,21,15,30,25,23,11,22,9,18);
@EXP64= (1,2,4,8,16,32,3,6,12,24,48,35,5,10,20,40,19,
	 38,15,30,60,59,53,41,17,34,7,14,28,56,51,37,
	 9,18,36,11,22,44,27,54,47,29,58,55,45,25,50,
	 39,13,26,52,43,21,42,23,46,31,62,63,61,57,49,33);	 

#@myarr=(1,1,1,1,1,1,1,1,1,1);
# print $GF16[17] ; 
# exit ;
sub tablecheck {
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	    printf "+ %2d" , &gfa($a,$b ) ;
	}
	print "\n" ; 
    }
	print "\n" ; 
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	    printf "* %2d" , &gfm($a,$b ) ;
	}
	print "\n" ; 
    }
	print "\n" ; 
}

sub vtablecheck {
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	    print " $a + $b = " , &gfa($a,$b ) , " : " ;
	}
	print "\n" ; 
    }
	print "\n" ; 
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	    print " $a * $b = " , &gfm($a,$b ), " : " ;
	}
	print "\n" ; 
    }
	print "\n" ; 
}

&set_fq();
if ( $exhaustive ) {
# Designed for larger q (>= 32) so we sample the space.
    $Max = 0.0 ; 
    $Sum=0.0;
    print "==== Sampling possibilities ====\n" ;
    for($samples=0;$samples<$Runs;$samples++){
	for($t=0;$t<$tr;$t++){
	    $myarr[$t]=int(1 + rand($Q));
	}
	if($init1){ $myarr[0]=1 };
      	&initialize(); $report="" ;
	for ( $t=0; $t<$tr ; $t++ ) {
	    $A = $myarr[$t];
	    &iterate();
	    $report .=  sprintf ("%d ", $A , $t, $H); 
	}
	&entropy() ; &scaleH();
	$Sum+=$H;
	if ( $verbose || ( $H > $Max * 0.9999 ) ) { 
	    print $report ; 	printf "Hscaled =%8g",  $H; 
	    print "\n" ;
	}
	if ( $H > $Max ) { $Max = $H ; @bestyet=@myarr }
    }
    printf "max was $Max , achieved by %s; average H=%8g \n", join(',',@bestyet),$Sum / $Runs; 
}
&say_channel() ;

##############################################################3
sub scaleH {
	if ( $q == 2 ) { $H =   2.0 * $H; }
	elsif ( $q == 8 ) { $H =  2.0 * $H / 3.0 ; }
	elsif ( $q == 16 ) { $H =  0.5 * $H ; }
	elsif ( $q == 32 ) { $H =  0.4 * $H ; }
	elsif ( $q == 64 ) { $H =  0.3333 * $H ; }
}
sub sayH {
	if ( $q == 2 ) { printf " 2 * "; $H =   2.0 * $H; }
	elsif ( $q == 8 ) { printf " 2/3 * "; $H =  2.0 * $H / 3.0 ; }
	elsif ( $q == 16 ) { printf " 1/2 * "; $H =  0.5 * $H ; }
	elsif ( $q == 32 ) { printf " 2/5 * "; $H =  0.4 * $H ; }
	elsif ( $q == 64 ) { printf " 1/3 * "; $H =  0.333 * $H ; }
	printf "H=%8g",  $H; 
	print "\n" ;
}
sub doit {
    &initialize();
    for ( $t=1; $t<=$tr ; $t++ ) {
	&iterate_chat() ;     }
}
sub doitA3 {
    &initialize();
    for ( $t=1; $t<=$tr ; $t++ ) {
	&iterate_chat() ; 
	$A ++ ; if ( $A > $Q ) { $A = $Q-1 ; }
    }
}
sub doitA8 {
    &initialize();
    for ( $t=1; $t<=$tr ; $t++ ) {
	&iterate_chat() ; 
	$A ++ ; $A *= 2 ; while ( $A > $Q ) { $A -= $Q-2 ; }
    }
}
sub doitA2 {
    &initialize();
    for ( $t=1; $t<=$tr ; $t++ ) {
	&iterate_chat() ; 
	$A *= 2 ; if ( $A > $Q ) { $A = 1 ; }
    }
 }
sub doitA {
    &initialize();
    for ( $t=1; $t<=$tr ; $t++ ) {
	&iterate_chat() ; 
	$A ++ ; if ( $A > $Q ) { $A = 1 ; }
    }
}
sub iterate_chat {
	&iterate();
	&entropy() ; 
	printf "A %d col %d H=%8g", $A , $t, $H; 
	if ( $q == 2 ) { printf " 2H=%8g",  2.0 * $H; }
	elsif ( $q == 8 ) { printf " H*2/3=%8g",  2.0 * $H / 3.0 ; }
	elsif ( $q == 16 ) { printf " H/2=%8g",  0.5 * $H ; }
	elsif ( $q == 32 ) { printf " 2H/5=%8g",  0.4 * $H ; }
	elsif ( $q == 64 ) { printf " H/3=%8g",  0.333 * $H ; }
	print "\n" ;}

sub iterate {
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	$sum = 0.0 ; 
	for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	    $sum += $p[&gfa(&gfm($b,$A),$a)] * $fq[$b] ;
	}
	$pnew[$a] = $sum ; 
    }
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	$p[$a] = $pnew[$a] ;
    }
}
sub gfa {
    local ( $a , $b ) = @_ ;
    return $a ^ $b ; 
#    if ( $a == 0 ) { return $b ; }
#    elsif ($a 
}

sub gfm {
    local ( $a , $b ) = @_ ;
    if ( $a == 0 || $b == 0 ) { return 0 ; }
    if ( $a == 1 ) { return $b ; }
    if ( $b == 1 ) { return $a ; }
    SWITCH: {
	if ( $q == 4 )  { return $EXP4[($LOG4[$a]+$LOG4[$b])%3]}
	if ( $q == 8 )  { return $EXP8[($LOG8[$a]+$LOG8[$b])%7]; }
	if ( $q == 16 )  { return $EXP16[($LOG16[$a]+$LOG16[$b])%15]; }
	if ( $q == 32 )  { return $EXP32[($LOG32[$a]+$LOG32[$b])%31]; }
	if ( $q == 64 )  { return $EXP64[($LOG64[$a]+$LOG64[$b])%63]; }
    }
    print "ERR" ; 
    return 0 ;
}

sub initialize {
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	$p[$a] = 0.0 ;
    }
    $p[0] = 1.0 ; 
}
sub entropy {
    $H = 0.0 ; 
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	$H -= &plgp( $p[$a] ) ;
    }
}
sub H2 {
    # binary entropy in natural units
    # usage ( $f ) 
    local ( $f ) = @_ ; 
    local ( $ans ) = - (  $f * &safelog ( $f ) + (1.0 - $f) * &safelog (1.0 - $f ) )  ;
    return $ans ; 
}

sub Binom {
    my ( $a , $x ) = @_;
    local ( $Bits, $i, $result, $y );

    $y=int(0.5+log($q)/log(2.0));
    $result=1;
    $Bits=0;
#   Count number of bits set in binary number `a'. 
    for($i=16;$i>0;$i--){
	if($a&0x1) {$Bits++;}
	$a>>=1;
    }
    for($i=0;$i<$Bits;$i++){
	$result*=$x;
    }
    for(;$i<$y;$i++){
	$result*=(1-$x);
    }
    return($result);
}

sub set_fq {
    if ( $qary ) { # uniform noise
	$qf = $f / $Q ; 
	for($i=1;$i<=$Q;$i++){ $fq[$i] = $qf ; }
	$fq[0] = 1 - $f ; 
    }
    else {
	for($i=0;$i<$q;$i++){
	    $fq[$i] = &Binom($i,$f);
	}
    }
    &say_channel() ;
}

sub say_channel {
    # report entropy of channel
    for ( $i = 0 ; $i <= $Q ; $i ++ ) { $p[$i] = $fq[$i] ; }
    &entropy();
    print "noise has entropy: " ;
    &sayH();
}
sub plgp {
    local ($x) = @_ ; 
    return $x * &safelog($x) / log(2.0) ; 
}

sub safelog {
    local ($f) = @_ ; 
    if ( $f < 0 ) { print STDERR "Can't take log of $f\n" ; }
    return ( $f > 0 ? log ( $f) : 0.0 ) ;
}