# Copyright 1999-2000 Mats Kindahl. All rights reserved.
# This program is free software; you can redistribute it and/or
# modify it under the same terms as Perl itself.
#
# E-mail: matkin@docs.uu.se
# 
# Program to compute the linear complexity of a bit sequence using the
# Berlekamp-Massay algorithm.

use strict;

sub print_poly {
    my $string = "1";
    foreach my $i (1..(@_-1)) {
	$string .= " + D^$i" if $_[$i] > 0;
    }
    return $string;
}

# The Berlekamp-Massay algorithm to compute the linear complexity
# profile of a binary sequence. The sequence is fed as argument to the
# function. The pseudo-code of the algorithm is as follows.
#
# C(D) := 1
# L := 0
# m := -1
# B(D) := 1
# for N := 0 to n-1 do
#    d := (S[n] + sum(i : 1 <= i <= L : C[i]*S[N-i])) mod 2
#    if d = 1 then
#       T(D) := C(D)
#       C(D) := C(D) + B(D) * D^(N-m)
#       if L <= N div 2 then
#          L := N + 1 - L
#          m := N
#          B(D) := T(D)
#       end if
#    end if
# end for 
# return [L,C(D)]

sub complexity_profile {
    my @C = (1);		# Current generating polynomial
    my $c_len = 0;		# Current complexity
    my @B = (1);		# Previous generating polynomial
    my $b_len = -1;		# Previous complexity

    foreach my $bit (0..(@_-1)) {
	my $discrepancy = $_[$bit];
        # Compute discrepancy
	foreach my $i (1..$c_len) { 
	    $discrepancy += $C[$i] * $_[$bit-$i];
	}
	$discrepancy = $discrepancy % 2;

	if ($discrepancy == 1) {
	    my @T = @C;
	    # C(D) = C(D) + B(D) ^ (N-m)
	    foreach my $i (0..(@B-1)) {
		$C[$i+($bit-$b_len)] += $B[$i];
		$C[$i+($bit-$b_len)] %= 2;
	    }
	    if ($c_len <= $bit / 2) {
		$c_len = $bit + 1 - $c_len;
		$b_len = $bit;
		@B = @T;
	    }
	}
    }
    return ($c_len, @C);
}

my $verbose = 0;

my $print_complexity = 1;
my $print_tap = 0;

my $line_mode = 1;

my $usage = <<"END_MESSAGE";
Usage: $0 <option>
--verbose         Be verbose (default off)
--unpack          Unpack input using 'B*'
--complexity=[01] Display complexity (default 1)
--tap=[01]        Display tap sequence (default 0)
--line=[01]       Read one line at a time instead of whole input (default 1).
END_MESSAGE

# Default unpack routine to use (no unpacking).
my $unpack = sub { return split("", $_[0]); };

# Read the options
OPTION: while (my $arg = shift(@ARGV)) {
    if ($arg !~ /^-/) {
	unshift(@ARGV, $arg);
	last OPTION;
    }

    if ($arg =~ /^--unpack/) {
	$unpack = sub { return split("", unpack('B*', $_[0])) };
    }
    elsif ($arg =~ /^--verbose/) {
	$verbose = 1;
    }
    elsif ($arg =~ /^--complexity=(\d+)/) {
	$print_complexity = $1;
    }
    elsif ($arg =~ /^--tap=(\d+)/) {
	$print_tap = $1;
    }
    elsif ($arg =~ /^--line=(\d+)/) {
	$line_mode = $1;
    }
    else {
	die $usage;
    }
}

sub do_computation (&) {
    if ($line_mode) {
	while (<>) {
	    $_[0]->($_);
	}
    } else {
	my @lines = <>;
	$_[0]->(join('',@lines));
    }
}

do_computation {
    my @sequence = $unpack->($_[0]);
    my ($complexity, @tap) = complexity_profile(@sequence);

    if ($verbose) {
	my $input = "Input (length " . length($_[0]) . ")";
	print $input, "\n", '-' x length($input), "\n", $_[0], "\n";

	my $tap_sequence = "Tap Sequence (linear complexity $complexity)";
	print "$tap_sequence\n";
	print '-' x length($tap_sequence), "\n";
	print join('', @tap), "\n";
    } else {
	print $complexity if $print_complexity;
	print "\t" if $print_complexity and $print_tap;
	print @tap if $print_tap;
	print "\n";
    }
}