/* @(#)b-m.c	1.3 (QUALCOMM) 09/02/02 */
/* Berlekamp-Massey algorithm */
/* see Menezes et al, P.201 */

#include <unistd.h>
#include <stdio.h>
#include <string.h>

typedef unsigned char uchar;
#define MAX 100000
int	verbose = 0;
int	polyform = 0;

uchar	inbuf[MAX / 8]; /* byte-binary input buffer */
int	n, N;		/* current / total number of input bits */
int	L;		/* current linear complexity */
int	m;
uchar	s[MAX], C[MAX], B[MAX], T[MAX];		/* polynomials */
uchar	d;		/* current discrepancy */

#define cp(x,y) memcpy(x,y,n+1);

void
pr(char *s, uchar *a)
{
    register int	i;
    int			first = 1;

    if (!polyform) {
	printf("%s", s);
	for (i = 0; i <= L; ++i) {
	    printf("%d", a[i]);
	}
    }
    else {
	for (i  = 0; i <= L; ++i) {
	    if (a[i]) {
		if (first)
		    first = 0;
		else
		    printf(" + ");
	    }
	    if (i == 0 && a[i])
		printf("1");
	    else if (i == 1 && a[i])
		printf("x");
	    else if (a[i])
		printf("x^%d", i);
	}
    }
    printf("\n");
}

int
main(int ac, char **av)
{
    int		i, j;
    char	*myname = av[0];
    int		warned = 0;

    while (ac >= 2 && av[1][0] == '-') {
	if (strcmp(av[1], "-v") == 0) {
	    verbose = 1;
	    ++av, --ac;
	}
	else if (strcmp(av[1], "-p") == 0) {
	    polyform = 1;
	    ++av, --ac;
	}
    }
    if (ac > 2) {
	fprintf(stderr, "Usage: %s [-v|-p] [bits]\n  (or binary on stdin)\n", myname);
	return 1;
    }
    else if (ac == 2) {
	N = strlen(av[1]);
	for (i = 0; i < N; ++i)
	    s[i] = av[1][i] == '0' ? 0 : 1;
    }
    else {
	N = fread(inbuf, 1, MAX/8, stdin);
	if (N <= 0) {
	    fprintf(stderr, "%s: no input\n", myname);
	    return 1;
	}
	/* expand inbuf to be one bit per byte */
	for (i = 0; i < N; ++i) {
	    for (j = 0; j < 8; ++j) {
		s[8*i + j] = (inbuf[i] >> (7 - j)) & 0x1;
	    }
	}
	N *= 8;
    }
    /* s[0..N] is the input sequence */

    /* initialise */
    n = 0;
    C[n] = B[n] = 1;
    L = 0;
    m = -1;

    /* accept a bit of the input sequence, adjusting as needed */
    while (n < N) {
	/* compute discrepancy d */
	d = s[n];
	for (i = 1; i <= L; ++i)
	    d ^= C[i] & s[n - i];
	/* recompute sequence if necessary */
	if (d) {
	    cp(T, C);
	    for (i = n-m, j = 0; i <= n+1; ++i, ++j)
		C[i] ^= B[j];
	    if (L <= n/2) {
		/* we have a jump in linear complexity -- check if unusual */
		if ((n/2 - L) > 6 && L > 2 && !warned) {
		    printf("Large jump (%d) in linear complexity detected\n",
			n - 2*L);
		    warned = 1;
		}
		L = n + 1 - L;
		m = n;
		cp(B, T);
	    }
	    if (verbose) {
		printf("n=%d s[%d]=%d d=%d L=%d m=%d\n", n, n-1, s[n-1], d, L, m);
		pr("T=", T);
		pr("C=", C);
		pr("B=", B);
	    }
	}
	++n;
    }
    /* count weight of polynomial -- don't include last bit */
    for (i = j = 0; i < L; ++i)
	j += C[i];
    printf("L=%d W=%d\n", L, j);
    pr("", C);
    if (L > (N/2-4)) {
	printf("Sequence is probably not linear");
	if (N == MAX)
	    printf(" up to MAX (%d) tested", MAX);
	printf("\n");
	return warned;
    }
    return 2;
}