#!/usr/local/bin/perl -w
#
#                                                  Free Software, (c) David MacKay August 2003
#                                                  Distributed under the GPL.
#                                                  http://www.gnu.org/copyleft/gpl.html
#
# hz.p program to find the entropy of z when it is given by A x = z 
# using GF(4)
#
# Written before Sep 10 1997
#
# reused 1998 when investigating the Richardson Urbanke fft method.
#
# @p used to hold the evolving probability distribution
# Now @P holds its evolving FFT.
#
# hzFFTe.p is the efficient version that only visits each possible
# row permutation once.
#
# Apologies, this code is very messy.
# If modifying it, I would recommend focussing on the way in which
# GF(64) is handled. The GF(64) code was added most recently and
# is most elegant.  For GF(8) I was "lazy" and coded the multiply
# and add tables by hand.  GF(16) also has explicit tables.
# In all cases the addition and multiplication tables assumed are
# written to STDOUT so they can be checked.

# usage:
# 
# hzFFTe.p tr=3 q=16 exhaustive=1 > tmpe16.3
# hzFFTe.p tr=4 q=16 exhaustive=1 > tmpe16.4
# hzFFTe.p tr=5 q=16 exhaustive=1 > tmpe16.5
# hzFFTe.p tr=6 q=16 exhaustive=1 > tmpe16.6
# hzFFTe.p tr=7 q=16 exhaustive=1 > tmpe16.7
#
# hzFFTe.p tr=3 q=64 exhaustive=1 > tmpe64.3
# hzFFTe.p tr=4 q=64 exhaustive=1 > tmpe64.4
# hzFFTe.p tr=5 q=64 exhaustive=1 > tmpe64.5
#
#  see code/GFQ/README


$Auniform = 0 ; 
$tr = 4 ; # number of elemnts to add. (row weight)
$f = 0.1 ; # noise of BSC
$q = 4 ;  # which field to do
$verbose = 0 ; 
$exhaustive = 0 ; 
$qary = 0 ; 
$doperm = 1 ; # whether to try some random perms, and how many

#  How to override with command line
eval "\$$1=\$2" while @ARGV && $ARGV[0]=~ /^(\w+)=(.*)/ && shift;

$GFQ = 1 ; # to start with, we are using GFQ. this gets switched off in doperm
$B =  $q == 4 ? 2 : ( $q == 8 ) ? 3 : ($q == 16) ? 4: ($q == 64) ? 6 : ($q==2)?1: 0  ;
if ( !$B ) { print "Hell baby I jus' don't know what this damn B is! $q\n" ;}

if ($B==6) {
    $Q = 64 ;
    $mask = "1000011" ;
    $mask = 2**6+2+1 ;
    &setupm;
    # creates look up table for multiply
}


$Q = $q-1 ; 
@GF16=(0,0,0,0, 0,0,0,0, 0 ,0 ,0 ,0 , 0 ,0 ,0 ,0 ,
0,1 ,2 ,3 , 4 ,5 ,6 ,7 , 8 ,9 ,10,11, 12,13,14,15,
0,2 ,4 ,6 , 8 ,10,12,14, 3 ,1 ,7 ,5 , 11,9 ,15,13,
0,3 ,6 ,5 , 12,15,10,9 , 11,8 ,13,14, 7 ,4 ,1 ,2 ,
0,4 ,8 ,12, 3 ,7 ,11,15, 6 ,2 ,14,10, 5 ,1 ,13,9 ,
0,5 ,10,15, 7 ,2 ,13,8 , 14,11,4 ,1 , 9 ,12,3 ,6 ,
0,6 ,12,10, 11,13,7 ,1 , 5 ,3 ,9 ,15, 14, 8, 2, 4,
0,7 ,14,9 , 15,8 ,1 ,6 , 13,10,3 ,4 , 2 , 5,12,11,
0,8 ,3 ,11, 6 ,14,5 ,13, 12,4 ,15,7 , 10, 2, 9, 1,
0,9 ,1 ,8 , 2 ,11,3 ,10, 4 ,13,5 ,12, 6 ,15, 7,14,
0,10,7 ,13, 14,4 ,9 ,3 , 15,5 ,8 ,2 , 1 ,11,6 ,12,
0,11,5 ,14, 10,1 ,15,4 , 7 ,12,2 ,9 , 13,6 ,8 ,3 ,
0,12,11,7 , 5 ,9 ,14,2 , 10,6 ,1 ,13, 15,3 ,4 ,8 ,
0,13,9 ,4 , 1 ,12,8 ,5 , 2 ,15,11,6 , 3 ,14,10,7 ,
0,14,15,1 , 13,3 ,2 ,12, 9 ,7 ,6 ,8 , 4 ,10,11,5 ,
0,15,13,2 , 9 ,6 ,4 ,11, 1 ,14,12,3 , 8 ,7 ,5 ,10)
    ;
# 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" ; 
}

if ( $verbose ) { &vtablecheck() ; } else { &tablecheck() ; }
&set_fq();
# find FFT of the channel distbn.
@FQ = &fft ( @fq ) ; 
if ( $verbose ) {
# check we can get it back...
    @fq2 = &fft ( @FQ ) ; 
    printf "%2s\t%5s\t%5s\t%5s\n" , "q", "f(q)" , "F(q)" , "f(q)" ; 
    for ( $qq = 0 ; $qq <= $Q ; $qq ++ ) {
	printf "%2d\t%5.2f\t%5.2f\t%5.2f\n" , $qq , $fq[$qq] , $FQ[$qq] , $fq2[$qq]/(1.0*$q) ;
    }
}
if ( $Auniform ) {
    for ($A=1;$A<=$Q;$A++) {
	print "==== all A are $A ====\n" ; 
	&doit() ; 
	if ( $qary ) { exit ; }
    }
}
if ( $q > 2 ) {
    for ( $k = $Q ; $k <= $Q ; $k ++ ) {
	$A = $k ; 
	print "==== A goes 1,2,..Q,1 e.g. 1,2,3,1,2,3 === " , $k , "\n" ; 
	&doitA() ; 
    }
    print "==== A goes Q-1,Q,Q-1,Q e.g. 3,2,3,2 ====\n" ; 
    $A = $Q ; 
    &doitA3() ; 
    if ( $q > 4 ) {
	print "==== A pseudorandom ====\n" ; 
	$A = $Q ; 
	&doitA8() ; 
	print "==== A goes 1,2,4,8.. ====\n" ; 
	$A = 1 ; 
	&doitA2() ; 
    }
}
sub install_new_As {
    local($t) = @_ ;
    $As[$t] ++ ;
    if ( $As[$t] <= $Q ) {
	# done 
    } elsif ($t==$tr-1) {
	$As[$t] = 0 ; # terminate exhaustive thing
    } else {
	$As[$t] = install_new_As($t+1) ;
    }
    return $As[$t] ; 
}
if ( $exhaustive ) {
    $Max = 0.0 ; 
    $done = 0 ; 
    $count=0 ; $Tot=0.0;
    # set up first one
    for ( $t=1; $t<=$tr ; $t++ ) {
	$As[$t] = 1 ;
    }
    $As[1] = 0 ; 
    $As[0] = "" ; 
    print "==== exhaustive enumerate\n" ;
    while ( install_new_As(1) ) {
	&initialize(); $report="" ; 
	for ( $t=1; $t<=$tr ; $t++ ) {
	    $A = $As[$t] ;  
	    &iterate();
	    $report .= sprintf "%d ", $A ; 
	}
	&fftentropy() ; &scaleH();
	if ( $verbose || ( $H > $Max * 0.9999 ) ) { 
	    print $report ; 	printf "Hscaled =%8g",  $H; 
	    print "\n" ;
	}
	if ( $H > $Max ) { $Max = $H ; $kmax = join(":", @As) ; }
	if ( $exhaustive != 1 ) { $done = 1 ; }
	$Tot += $H ; $count ++ ; 
    } $Mean = $Tot / $count ; 
    print "max was $Max , achieved by $kmax| \n" ; 
    print "mean of $count was $Mean (above list gives record-breakers only)\n";
}
if ( $doperm ) {
    $GFQ = 0 ; 
    $Max = 0.0 ; 
    $count=0 ; $Tot=0.0;
    print "==== permutations instead of GF(q)\n" ; 
    for ( $k = 1 ; 
	 $k <= $doperm ;
	 $k ++ ) {
	&initialize(); $report="$k| " ; 
	for ( $t=1; $t<=$tr ; $t++ ) {
	    @perm = &permutation ( $Q , 0 ) ;
	    print join ( ':', @perm) . "\n" if ( $verbose >=3 )  ; 
	    &iterate();
	    $report .= join(':',@perm) ." " ; 
	}
	&fftentropy() ; &scaleH();
	if ( $verbose || ( $H > $Max * 0.9999 ) ) { 
	    # only reports pretty good answers
	    print $report ; 	printf "Hscaled =%8g",  $H; 
	    print "\n" ;
	}
	if ( $H > $Max ) { $Max = $H ; $kmax = $k ; }
	$Tot += $H ; $count ++ ; 
    } $Mean = $Tot / $count ; 
    print "max was $Max , achieved by $kmax| \n" ; 
    print "mean was $Mean\n";
}
&say_channel() ;

##############################################################3
sub scaleH { # express in terms comparable to GF(4)
	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 == 64 ) { $H =  2.0 * $H / 6.0 ; }
}
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 == 64 ) { printf " 2/6 * "; $H =  2.0 * $H / 6.0 ; }
	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();
	&fftentropy() ; 
	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 == 64 ) { printf " H*2/6=%8g",  2.0 * $H / 6.0 ; }
	print "\n" ;}

# usage @perm = &permutation ( $Q , 0 ) ;
sub permutation { # returns a permn of 0..Q offset by offset eg 0/1
    local ( $Q , $offset ) = @_ ;
    for ( local($q) = 0 ; $q <= $Q ; $q ++ ) { $qs[$q] = $q + $offset ; }
    for ( local($n) = 0 ; $n <= $Q ; $n ++ ) {
#	$from = $n ;
	$q = int ( rand($#qs+1) ) ;
	$to[$n] = splice ( @qs , $q , 1  ) ;
	print "to[$n] = $to[$n]\n" if ( $verbose >=3 ) ; 
    }
    print "Here is a perm $Q:\n"  if ( $verbose >=3 ) ; 
    print join ( ':' ,  @to  ) , "\n" if ( $verbose >=3 )  ; 
    return ( @to ) ; 
}

sub iterate { # maps p and fq to p using A 
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	if ( $GFQ ) {
	    $thisfq[&gfm($a,$A)] = $fq[$a] ; # this is the only place gfm is used
	} else {
	    $thisfq[$perm[$a]] = $fq[$a] ; 
	}
    } 
    @Thisfq = &fft( @thisfq ) ; # maybe this is wasteful; there is 
# probably some simple relatinoship between all the different FTs
# for the different values of $A. So we could 
# just fft @fq once only and then permute? Not quite that simple, there
# are sign changes too.
    if ( $verbose >= 2) {
	printf "%2s\t%5s\t%5s\n" , "q", "f(q)" , "F(q)"  ; 
	for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	    printf ("%2d\t%5.2f\t%5.2f\n" , $b , $thisfq[$b] , $Thisfq[$b] ) ;
	}
    }
    for ( $b = 0 ; $b <= $Q ; $b ++ ) {
	$P[$b] = $P[$b] * $Thisfq[$b] ;
    }
}
sub gfa {
    local ( $a , $b ) = @_ ;
    return $a ^ $b ; 
}
# example FFT: from (0.9,0.1)^4
#  q       f(q)    F(q)
#  0       0.66    0.41
#  1       0.07   -0.51
#  2       0.07   -0.51
#  3       0.01    0.64
#  4       0.07   -0.51
#  5       0.01    0.64
#  6       0.01    0.64
#  7       0.00   -0.80
#  8       0.07   -0.51
#  9       0.01    0.64
# 10       0.01    0.64
# 11       0.00   -0.80
# 12       0.01    0.64
# 13       0.00   -0.80
# 14       0.00   -0.80
# 15       0.00    1.00
sub fft { # receives a list argument, returns a list
#
# Note that this should divide by sqrt($q) if it is to be
# a self-inverse transform.
# I prefer to do a division by $q by hand at a convenient
# moment.
#
    local ( @p ) = @_  ;
    local ( $Q2 ) = $q/2 - 1 ;
    local ( $b )  = 1 ; local ( $factor ) = 1 ; 
    for (  ; $b <= $B ; $b ++ ) {
	print "b = $b\n" if ($verbose >= 3 ) ; 
	for ( $rest = 0 ; $rest <= $Q2 ; $rest ++ ) {
	    $restH = int($rest / $factor ) ; # integer
	    $restL = $rest % $factor ; 
	    print "$restH:$restL\n" if ($verbose >= 3 ) ; 
	    $rest0 = $restH * $factor * 2 + $restL ; 
	    $rest1 = $rest0 + $factor ; 
	    print "FFT munging $rest0,$rest1: $p[$rest0],$p[$rest1] to " if ( $verbose >= 3) ;
	    local ($sum) = $p[$rest0] + $p[$rest1] ;
	    $p[$rest0] =  $p[$rest1] - $p[$rest0] ;
	    $p[$rest1] =  $sum ; 
	    print "$p[$rest0],$p[$rest1]\n" if ( $verbose >= 3 ) ;
	}
	$factor = $factor * 2 ; 
    }
    return @p ; 
}

sub gfm {
    local ( $a , $b ) = @_ ;
    if ( $q == 64 ) { return $gf64{$a,$b} ; }
    elsif ( $q == 16 ) { return $GF16[ ( $a ) * 16 + ( $b ) ] ; }
    if ( $a == 0 || $b == 0 ) { return 0 ; }
    elsif ( $a == 1 ) { return $b ; }
    elsif ( $b == 1 ) { return $a ; }
    elsif ( $q == 4 ) {
	if ( $a == 2 ) {
	    if ( $b == 2 ) { return 3 ; }
	    elsif ( $b == 3 ) { return 1 ; }
	    else { print "ERR" ; }
	}
	elsif ( $a == 3 ) {
	    if ( $b == 2 ) { return 1 ; }
	    elsif ( $b == 3 ) { return 2 ; }
	    else { print "ERR" ; }
	}
	else { print "ERR" ; }
    }
    elsif ( $q == 8 ) {
	if ( $a == 2 ) {
	    if ( $b == 2 )    { return 4 ; }
	    elsif ( $b == 3 ) { return 6 ; }
	    elsif ( $b == 4 ) { return 3 ; }
	    elsif ( $b == 5 ) { return 1 ; }
	    elsif ( $b == 6 ) { return 7 ; }
	    elsif ( $b == 7 ) { return 5 ; }
	    else { print "ERR" ; }
	}
	elsif ( $a == 3 ) {
	    if ( $b == 2 )    { return  6; }
	    elsif ( $b == 3 ) { return  5; }
	    elsif ( $b == 4 ) { return  7; }
	    elsif ( $b == 5 ) { return  4; }
	    elsif ( $b == 6 ) { return  1; }
	    elsif ( $b == 7 ) { return  2; }
	    else { print "ERR" ; }
	}
	elsif ( $a == 4 ) {
	    if ( $b == 2 )    { return  3; }
	    elsif ( $b == 3 ) { return  7; }
	    elsif ( $b == 4 ) { return  6; }
	    elsif ( $b == 5 ) { return  2; }
	    elsif ( $b == 6 ) { return  5; }
	    elsif ( $b == 7 ) { return  1; }
	    else { print "ERR" ; }
	}
	elsif ( $a == 5 ) {
	    if ( $b == 2 )    { return 1 ; }
	    elsif ( $b == 3 ) { return 4 ; }
	    elsif ( $b == 4 ) { return 2 ; }
	    elsif ( $b == 5 ) { return 7 ; }
	    elsif ( $b == 6 ) { return 3 ; }
	    elsif ( $b == 7 ) { return 6 ; }
	    else { print "ERR" ; }
	}
	elsif ( $a == 6 ) {
	    if ( $b == 2 )    { return 7 ; }
	    elsif ( $b == 3 ) { return 1 ; }
	    elsif ( $b == 4 ) { return 5 ; }
	    elsif ( $b == 5 ) { return 3 ; }
	    elsif ( $b == 6 ) { return 2 ; }
	    elsif ( $b == 7 ) { return 4 ; }
	    else { print "ERR" ; }
	}
	elsif ( $a == 7 ) {
	    if ( $b == 2 )    { return 5 ; }
	    elsif ( $b == 3 ) { return 2 ; }
	    elsif ( $b == 4 ) { return 1 ; }
	    elsif ( $b == 5 ) { return 6 ; }
	    elsif ( $b == 6 ) { return 4 ; }
	    elsif ( $b == 7 ) { return 3 ; }
	    else { print "ERR" ; }
	}
	else { print "ERR" ; }
    }
    else { print "ERR" ; }
    return 0 ;
}

sub initialize {
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	$p[$a] = 0.0 ;
    }
    $p[0] = 1.0 ; 
    @P = &fft ( @p ) ; 
}
sub fftentropy {
    @p = fft( @P ) ; 
    for ( $a = 0 ; $a <= $Q ; $a ++ ) {
	$p[$a] = $p[$a] / $q ;
    }
    &entropy();
}
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 set_fq {
    if ( $qary ) { # uniform noise
	$qf = $f / $Q ; 
	for($i=1;$i<=$Q;$i++){ $fq[$i] = $qf ; }
	$fq[0] = 1 - $f ; 
    }
    else {
	if ( $q == 8 ) {
	    $fq[0] = (1-$f)*(1-$f)*(1-$f) ; 
	    $fq[1] = $f*(1-$f)*(1-$f) ; 
	    $fq[2] = (1-$f)*$f*(1-$f) ; 
	    $fq[3] = ($f*$f)*(1-$f) ; 
	    $fq[4] = (1-$f)*(1-$f)*($f) ; 
	    $fq[5] = $f*(1-$f)*($f) ; 
	    $fq[6] = (1-$f)*$f*($f) ; 
	    $fq[7] = ($f*$f)*($f) ; 
	} elsif ( $q == 16 ) {
	    $fq[0] = (1-$f)*(1-$f)*(1-$f)*(1-$f) ; 
	    $fq[1] = $f*(1-$f)*(1-$f)*(1-$f) ; 
	    $fq[2] = (1-$f)*$f*(1-$f)*(1-$f) ; 
	    $fq[3] = ($f*$f)*(1-$f)*(1-$f) ; 
	    $fq[4] = (1-$f)*(1-$f)*($f)*(1-$f) ; 
	    $fq[5] = $f*(1-$f)*($f)*(1-$f) ; 
	    $fq[6] = (1-$f)*$f*($f)*(1-$f) ; 
	    $fq[7] = ($f*$f)*($f)*(1-$f) ; 
	    $fq[8] = (1-$f)*(1-$f)*(1-$f)*($f) ; 
	    $fq[9] = $f*(1-$f)*(1-$f)*($f) ; 
	    $fq[10] = (1-$f)*$f*(1-$f)*($f) ; 
	    $fq[11] = ($f*$f)*(1-$f)*($f) ; 
	    $fq[12] = (1-$f)*(1-$f)*($f)*($f) ; 
	    $fq[13] = $f*(1-$f)*($f)*($f) ; 
	    $fq[14] = (1-$f)*$f*($f)*($f) ; 
	    $fq[15] = ($f*$f)*($f)*($f) ; 
	} elsif ( $q == 64 ) {
	    for ( $i = 0 ;$i <= $q-1 ; $i ++ ) {
		$fq[$i] = 1 ; $ii = $i ;
		for ( $b = 0 ; $b < $B ; $b ++ ) {
		    if ( $ii % 2 ) { $fq[$i] *= $f ; }
		    else { $fq[$i] *= (1-$f) ; }
		    $ii /= 2 ; 
		}
	    }
	} elsif ( $q == 4 ) {
	    $fq[0] = (1-$f)*(1-$f) ; 
	    $fq[1] = $f*(1-$f) ; 
	    $fq[2] = (1-$f)*$f ; 
	    $fq[3] = ($f*$f) ; 
	} elsif ( $q == 2 ) {
	    $fq[0] = 1-$f ; 
	    $fq[1] = $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 ) ;
}

sub setupm {
    $i=0;  # enter the look up table
    for ( $j = 0 ; $j < $Q ; $j ++ ) {
	$gf64{$j,0}=0;
	$gf64{0,$j}=0;
    }
    for ( $j = 1 ; $j < $Q ; $j ++ ) {
	for ($b=0;$b<$B;$b++) {
	    # create j shifted to the left by 0...B-1
	    $s[$b] = $j << $b ;
	}
	for ( $i = 1 ; $i < $Q ; $i ++ ) {
	    $m=0 ; $ii = $i ;
	    for ($b=0;$b<$B;$b++) {
		if ( $ii % 2 ) {
		    $m ^= $s[$b] ;
		}
		$ii = $ii / 2 ; # shift right and lose the bit.
	    }
	    # now subtract out all the overflowing stuff.
	    for ($b=$B-1;$b>=0;$b--) {
		if ( ( $m & ( 1 << ($b+$B) ) ) > 0 ) {
		    $m ^= ($mask << $b) ;
		}
		# assert...
		if ( ( ( $m & ( 1 << ($b+$B) ) ) > 0 ) ) {
		    print "something wrong!\n" ; 
		}
	    }
	    $gf64{$i,$j} = $m ;
	}
    }
    if($verbose) {
	print "\n" ; 
	for ( $j = 1 ; $j < $Q ; $j ++ ) {
	    for ( $i = 1 ; $i <= $j ; $i ++ ) {
		printf "%2d " , 	    $gf64{$i,$j}  ;
		if ( !($gf64{$i,$j} == $gf64{$j,$i} ) ) {
		    print "error!\n" ; 
		    printf "%2d,%2d " ,$gf64{$i,$j},$gf64{$j,$i} ;
		}
	    }
	    print "\n" ; 
	}
    }
    print "\n" ; 
}