3507.sagews

w=6; # Width of data to check
k=3; # Number of multiples of w to store as state (greater numbers increases safety)
n=w*k; # Number of bits of state

# Find some maximal length matrix for the underlying generator,
# so that if we get a sequence of zeros (for example), then the
# system keeps going
def make_A(t):
    while true:
        M=matrix(GF(2),n,n);
        for i in range(0,n):
            M[i,(i+w-1)%n]=1;
        for i in range(0,w):
            for j in range(1,t):
                M[i,randint(0,n-1)]=1;
        if not M.is_invertible():
            continue;
        if M.multiplicative_order()==2^n-1:
            return M;

# Find another matrix such that [A^(k-1)B ; A^(k-2)B;...;A*B;B] is invertible
def make_B(A,t=1):
    while true:
        B=matrix(GF(2),n,w);
        p=Permutations(w).random_element();
        for i in range(0,w):
            B[i,p[i]-1]=1;
        T=matrix(GF(2),n,0);
        F=identity_matrix(GF(2),n);
        for i in range(0,k):
            T=T.augment(F*B);
            F=F*A;
        if not T.is_invertible():
            return B;

def dump_taps(A,B):
    for i in range(0,n):
        print "x[%d] = 0 "%i,
        for j in range(0,n):
            if A[i][j]:
                print "^x[%d]"%j,
        for j in range(0,w):
            if B[i][j]:
                print "^s[%d]"%j,
        print ";";

A=make_A(3); print A.str()

[0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0] [0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0] [0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0] [0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0] [0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0] [0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1] [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0]

B=make_B(A,1); B

[0 0 1 0 0 0] [0 0 0 0 1 0] [0 0 0 1 0 0] [0 0 0 0 0 1] [1 0 0 0 0 0] [0 1 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 0 0 0]

dump_taps(A,B)

x[0] = 0 ^x[5] ^x[13] ^x[16] ^s[2] ; x[1] = 0 ^x[6] ^x[12] ^s[4] ; x[2] = 0 ^x[7] ^x[13] ^x[16] ^s[3] ; x[3] = 0 ^x[8] ^x[9] ^x[12] ^s[5] ; x[4] = 0 ^x[5] ^x[6] ^x[9] ^s[0] ; x[5] = 0 ^x[3] ^x[10] ^s[1] ; x[6] = 0 ^x[11] ; x[7] = 0 ^x[12] ; x[8] = 0 ^x[13] ; x[9] = 0 ^x[14] ; x[10] = 0 ^x[15] ; x[11] = 0 ^x[16] ; x[12] = 0 ^x[17] ; x[13] = 0 ^x[0] ; x[14] = 0 ^x[1] ; x[15] = 0 ^x[2] ; x[16] = 0 ^x[3] ; x[17] = 0 ^x[4] ;

def make_data(count):
    return random_matrix(GF(2),w,count);

def corrupt_data(data,n):
    M=matrix(GF(2),data.nrows(),data.ncols());
    G=copy(data);
    for i in range(0,n):
        while true:
            r=randint(0,data.nrows()-1);
            c=randint(0,data.ncols()-1);
            if M[r,c]==0:
                break;
        M[r,c]=1;
        G[r,c]=not G[r,c];
    return G;
    
def count_errors(data,got):
    return len((data-got).nonzero_positions());

def run_check(A,B,data):
    s=vector(GF(2),n);
    for i in range(0,n):
        s[i]=1;
    for i in range(0,data.ncols()):
        s=A*s+B*data.column(i);
    return s;

def calc_states(A,B,data):
    s=vector(GF(2),n);
    for i in range(0,n):
        s[i]=1;
    res=matrix(GF(2),n,0);
    res=res.augment(s);
    for j in range(0,data.ncols()):
        s=A*s+B*data.column(j);
        res=res.augment(s);
    return res;
    
def run_trial(A,B,count,nerrs):
    D=make_data(count);
    G=corrupt_data(D,nerrs);
    rs=run_check(A,B,D);
    gs=run_check(A,B,G);
    return rs==gs;

def plot_fails(A,B,count,reps):
    points=[];
    for nerrs in range(0,2*n):
        acc=0;
        for rep in range(0,reps):
            if run_trial(A,B,count,nerrs):
                acc=acc+1;
        points.append([nerrs,acc/reps]);
    return list_plot(points);

plot_fails(A,B,1000,n)