# Linear-Feedback Shift Register (LFSR)

In an LFSR, the output from a standard shift register is fed back into its
input causing an endless cycle. The feedback bit is the result of a linear
combination of the shift register content and the feedback coefficients.



The foolowing are the equation ruling the LFSR block scheme:
$$
\begin{align}
    b[t]       = &\; s_0[t] \\
    s_j[t]     = &\; s_{j+1}[t-1] \\
    s_{m-1}[t] = &\; \oplus_{j=0}^{m-1} p_{m-j} \otimes s_{j}[t-1] \\
\end{align}
$$

**Example**:
- LFSR length: 3
- feeedback polynomial: $x^3 + x + 1$
- initial state: 0b111
- LFSR cycle:
|     state | output | feedback |
|:---------:|:------:|:--------:|
| 0b111 (7) |    1   |     0    |
| 0b011 (3) |    1   |     1    |
| 0b101 (5) |    1   |     0    |
| 0b010 (2) |    0   |     0    |
| 0b001 (1) |    1   |     1    |
| 0b100 (4) |    0   |     1    |
| 0b110 (6) |    0   |     1    |
| 0b111 (7) |    1   |     0    |


```python
import math
from operator import xor
from functools import reduce
from itertools import compress
from itertools import islice
```

## LFSR generator

As first step, we implement an LFSR as a generator.

**Inputs**:
 - **Feedback Polynomial**: `list` of `int` representing the degrees of the non-zeros coefficients. Example: [12, 6, 4, 1, 0] represents the polynomial $x^{12}+x^6+x^4+x+1$.
 - **shift-register initial state** (optional, default all bits to 1): `int` or `list` of bits representing the LFSR initial state. Example: 0xA65 for 0b101001100101.
 
**Yield**:
- **Output bit**: `bool` representing the LFSR output bit

**Template**:
```python
def lfsr_generator(poly, state=None):
    ''' generator docstring '''
    
    # check inputs validity
    # define variables storing the internal state
    
    while True:
        # LFSR iteration:
        # - update state
        # - compute output
        yield output
```

### implentation with lists

In this first implementation both the shift register `state` and the feedback polynomial `poly` are `list` of `bool`. This is the most straightforward choice as it directly maps the LFSR block scheme.

| variable/<br>operation | definition | implementation |
|:-:|:-:|:--|
| feeedback polynomial | $$x^3 + x + 1$$ | `poly = [True, True, False, True]` |
| initial state | 0b111 | `state = [True, True, True]` |
| state update | $$s_j[t] = s_{j+1}[t-1]$$ | `state = state[1:] + [feedback]` |
| output bit | $$s_0[t]$$ | `output = state[0]` |
| feedback bit | $$\oplus_{j=0}^{m-1} p_{m-j} \otimes s_{j}[t]$$ | `feedback = reduce(xor, compress(state[::-1], poly[1:]))` |

For debugging purposes, a `verbose` flag is added as argument to enable the print of the internal state at each iteration.


```python
def lfsr_generator(poly, state=None, verbose=False):
    '''
    Generator implementing a Linear Feedback Shift Register (LFSR)
    
    Parameters
    ----------
    poly: list of int,
        feedback polynomial expressed as list of integers corresponding to
        the degrees of the non-zero coefficients.
    state: int, optional (default=None),
        shift register initial state. If None, all bits of state are set to 1.
    verbose: bool, optional (default=False),
        If True, the internal state is printed at each iteration
    
    Return
    ------
    bool, LFSR output bit
    '''
    
    # ==== check inputs ====
    # ==> verbose
    if verbose:
        # define a function to print the lfsr internal state
        def print_lfsr(state, output, feedback):
            state_bits = ''.join(str(int(s)) for s in state)[::-1]
            state_int = int(state_bits, 2)
            print(f'{state_bits} ({state_int})  {int(output)}  {int(feedback)}')
    
    # ==> poly
    length = max(poly) # LFSR length as the max degree of the polynomial
    poly = [i in poly for i in range(length+1)] # turn poly into a list of bool
    
    # ==> state
    if state is None: # default value for state is all ones (True)
        state = [True for _ in range(length)]
    else: # convert integer into list of bool
        state = [bool(int(s)) for s in (f'{state:0{length}b}')[::-1]]
    
    # ==== initial state ====
    # compute output bit and feedback bit
    output = state[0]
    feedback = reduce(xor, compress(state[::-1], poly[1:]))
    if verbose: # print initial state
        print(' state   b fb')
        print_lfsr(state, output, feedback)
    
    # ==== infinite loop ====
    while True: 
        state = state[1:] + [feedback]                          # update state 
        output = state[0]                                       # output bit
        feedback = reduce(xor, compress(state[::-1], poly[1:])) # feedback bit 
        if verbose: # print current state
            print_lfsr(state, output, feedback)
        
        yield output # return current output
```

Here we check the generator functioning, by replication the example showed above where the feedback polynomial is $x^3+x+1$ and the initial state is `0x7`=`0b111`


```python
poly = [3, 1, 0] # feedback polynomial x^3 + x + 1
state = 0x7      # initial state
niter = 7        # number of iterations

# lfsr definition
lfsr = lfsr_generator(poly, state, verbose=True) 

# just iter over the LFSR generator
for b in islice(lfsr, niter):
    # do nothing, just let the generator print its internal state
    pass
```

     state   b fb
    111 (7)  1  0
    011 (3)  1  1
    101 (5)  1  0
    010 (2)  0  0
    001 (1)  1  1
    100 (4)  0  1
    110 (6)  0  1
    111 (7)  1  0


### implentation with integers

In this implementation both the shift register `state` and the feedback polynomial `poly` are `int`. With this choice, bit-wise logical operation, as well as bit-shift, are easy to perform, while XOR of multiple bits or reversing the bit order are less straightforward.

For efficiency reasons, it is choosen to store the bits of the state in reverse order. At each iteration, the computation of the feedback bit requires a bit-wise AND between the current state and the polynomial feedback. Since these quantities are defined with inverted index, it is convinient to keep stored the bits of either state or poly in reverse order. To be compliant with common implementation in other programming languages (like C/C++) we choose to store the bits of the state in reverse order.

| variable/<br>operation | &nbsp; &nbsp; &nbsp; &nbsp; definition &nbsp; &nbsp; &nbsp; &nbsp; | implementation |
|:-:|:-:|:--|
| feeedback <br> polynomial | &nbsp; $x^3 + x + 1 $ | `poly = 0b1011`, `length = 3` |
| initial state | 0b111 | `state = 0b111` |
| state update | &nbsp; $s_j[t] = s_{j+1}[t-1] $ | `statemask = (1 << length) - 1` <br> `state = ((state << 1) \| feedback) & statemask` |
| output bit | &nbsp; $s_0[t] $ | `outmask = 1 << (length-1)` <br> `output = bool(state & outmask)` |
| feedback bit | &nbsp; $\oplus_{j=0}^{m-1} p_{m-j} \otimes s_{j}[t] $ &nbsp; | `feedback = parity(state & poly)` |

Let us define a function `reverse` to reverse the bit order of an integer


```python
def reverse(x, nbit=None):
    ''' reverse the bit order of the input `x` (int) represented with `nbit` 
    (int) bits (e.g., nbit: 5, x: 13=0b01101 -> output: 22=0b10110) '''
    if nbit is None:
        nbit = math.ceil(math.log2(x))
    return int(f'{x:0{nbit}b}'[::-1][:nbit], 2)
```


```python
x, nbit = 0b1101011000101010, 16
y = reverse(x, nbit)
print(f'{x:0{nbit}b} -> {y:0{nbit}b}')
```

    1101011000101010 -> 0101010001101011


Let us also define a function `parity` that computes the parity bit for a given integer.


```python
def parity(x):
    ''' compute the parity bit of an integer `x` (int) '''
    return bool(sum(int(b) for b in f'{x:b}') % 2)
```


```python
x = 0b1101011010101010
y = parity(x)

print(f'{x:0{nbit}b} -> parity bit: {y:b}')

```

    1101011010101010 -> parity bit: 1


Now we can define the generator that implements the LFSR.


```python
def lfsr_generator(poly, state=None, verbose=False):
    '''
    Generator implementing a Linear Feedback Shift Register (LFSR)
    
    Parameters
    ----------
    poly: list of int,
        feedback polynomial expressed as list of integers corresponding to
        the degrees of the non-zero coefficients.
    state: int, optional (default=None),
        shift register initial state. If None, all bits of state are set to 1.
    verbose: bool, optional (default=False),
        If True, the internal state is printed at each iteration
    
    Return
    ------
    bool, LFSR output bit
    '''
    
    # ==== check inputs ====
    # ==> verbose
    if verbose:
        def print_lfsr(state, output, feedback, length):            
            print(f'{state:0{length}b} ({state})',
                  f'{int(output)}  {int(feedback)}')
    
    # ==> poly
    length = max(poly) # LFSR length as the max degree of the polynomial 
    poly = sum([2**p for p in poly]) >> 1  # ignoring p0 (always 1)
    
    # ==> state
    outmask = 1 << (length-1)      # mask to select the output bit from state
    statemask = (1 << length) - 1  # mask for state
    if state is None: # default value for state is all ones
        state = statemask
    state = reverse(state, length) # reverse the bit order of state
    
    # ==== initial state ====
    # compute output bit and feedback bit
    output = bool(state & outmask) # get output from current state
    feedback = parity(state & poly)  # compute feedback bit as the parity bit
    if verbose: # print initial state
        print(' state   b fb')
        print_lfsr(reverse(state, length), output, feedback, length)
    
    # ==== infinite loop ====
    while True:
        state = ((state << 1) | feedback) & statemask # update state
        output = bool(state & outmask)                # get output
        feedback = parity(state & poly)               # update feedback bit 
        if verbose: # print current state
            _state = reverse(state, length)
            print_lfsr(_state, output, feedback, length)
        
        yield output
```

As we did above, we check the implementation by replicating the reference example.


```python
poly = [3, 1, 0]
state = 0x7
niter = 7

lfsr = lfsr_generator(poly, state, verbose=True)

for b in islice(lfsr, niter):
    pass
```

     state   b fb
    111 (7) 1  0
    011 (3) 1  1
    101 (5) 1  0
    010 (2) 0  0
    001 (1) 1  1
    100 (4) 0  1
    110 (6) 0  1
    111 (7) 1  0


Now, we can use the LFSR to generate a "long" (not really) sequence of bits.


```python
nbits = 1000
poly = [5, 2, 0]
lfsr = lfsr_generator(poly)

# iter over LFSR and store its output in a list
bits = [b for b in islice(lfsr, nbits)]

# print the list of bool as string of 0s and 1s
print(''.join([str(int(b)) for b in bits])) 
```

    1111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011010010000101011101100011111001101001000010101110110001111100110100100001010111011000111110011


We can notice that the 31-bits sequence `1111001101001000010101110110001` is periodically repeated. This is compliant with an LFSR with length 5 and a primitive polynomial with degree 5 in the feedback.

## LFSR Iterator

Implementing the LFSR as generator does not allow the user to access to the internal state and variables. Indeed, we needed to code a internal print that was activate by the `verbose` flag for debug.

To solve this issue, we can implement the LFSR as an **iterator**, that is a class with the `__iter__` and `__next__` methods that allow to iter over it. Such solution allow the user to employ the LFSR with the same scheme used for the generator bu with the capability to access to the internal state and variables.

**Inputs**:
 - **Feedback Polynomial**: `list` of `int` representing the degrees of the non-zeros coefficients. Example: [12, 6, 4, 1, 0] represents the polynomial $x^{12}+x^6+x^4+x+1$.
 - **shift-register initial state** (optional, default all bits to 1): `int` or `list` of bits representing the LFSR initial state. Example: 0xA65 for 0b101001100101.
 
**Attributes**:
- `poly`: (`list` of `int`) list of the polynomial coefficients;
- `length`: (`int`) polynomial degree and length of the shift register;
- `state`: (`int`) LFSR state;
- `output`: (`bool`) output bit;
- `feedback`: (`bool`) feedback bit.
 
**Methods**:
- `__init__`: class constructor;
- `__iter__`: necessary to be an iterable;
- `__next__`: update LFSR state and returns output bit;
- `cycle`: returns a list of bool representing the full LFSR cycle ;
- `run_steps`: execute N LFSR steps and returns the corresponding output as list of bool (N is a input parameter, default N=1);
- `__str__`: return a string describing the LFSR class instance.


**Template**:
```python
class LFSR(object):
    ''' class docstring '''

    def __init__(self, poly, state=None):
        ''' constructor docstring '''
        ...
        self.poly = ... 
        self.length = ... 
        self.state = ... 
        self.output = ... 
        self.feedback = ...

    def __iter__(self): 
        return self

    def __next__(self): 
        ''' next docstring ''' 
        ... 
        return self.output

    def run_steps(self, N=1): 
        ''' run_steps docstring ''' 
        ... 
        return list_of_bool

    def cycle(self, state=None): 
        ''' cycle docstring ''' 
        ... 
        return list_of_bool
```

We first define two functions to convert integers into, what we will call it, their *sparse* representation, which consist in the list of indexes corresponding to 1s in the integer binary representation. For example the integer 138 (`0b10001010`) is turned into the list `[1, 3, 7]`.

This function is employed to easily pass from the integer representation of the LFSR feedback polynomial (used during computations) to the representation as list of the degrees of the non-zeor coefficients (used as inputs and outputs).

We call it sparse as it reminds how sparse matrices are stored, i.e., as list of triplets comprising the row and column indexes and the corresponding value.


```python
def int_to_sparse(integer):
    ''' transform an integer into the list of indexes corresponding to 1
    in the binary representation (sparse representation)'''
    sparse = [i for i, b in enumerate(f'{integer:b}'[::-1]) if bool(int(b))]
    return sparse

def sparse_to_int(sparse):
    ''' transform a list of indexes (sparse representation) in an integer whose
    binary representation has 1s at the positions corresponding the indexes and 
    0 elsewhere '''
    integer = sum([2**index for index in sparse])
    return integer
```


```python
integer = 0b10001010

sparse = int_to_sparse(integer)
print(f'int to sparse: {integer} ({bin(integer)}) -> {sparse}')

integer = sparse_to_int(sparse)
print(f'sparse to int: {sparse} -> {integer} ({bin(integer)})')
```

    int to sparse: 138 (0b10001010) -> [1, 3, 7]
    sparse to int: [1, 3, 7] -> 138 (0b10001010)


Here is the LFSR implemented as Iterator. We choose to implement both state and polynomial feedback as integers.

Note that:
- in Python, for attributes and methods there is not the concept of **private** or **public** . It is all public and it is up to the user make a good use of the class. However, there is the convention (supported also by IDE) that attributes starting with an underscore `_` are to be considered private:
```python
self.state # this is public
self._state # this is considered private (actually, it is public)
```


- Python meets the idea of [Uniform Access Principle](https://en.wikipedia.org/wiki/Uniform_access_principle) that states *all services offered by a module should be available through a uniform notation, which does not betray whether they are implemented through storage or through computation*. 
    - This means that `lfsr.state` and `lfsr.state = 0` are much preferred than `lfsr.get_state()` and `lfsr.set_state(0)`. 
    - If you need to wrap the accesses to the attributes inside methods (to deny setting, to check the value before setting, to return the value in a different format with respect to how it is stored internally) you can use the [decorator](https://wiki.python.org/moin/PythonDecorators) `@property` that allows you to hide a method `lfsr.set_state(0)` behind the expression `lfsr.state = 0`.
    

These features help in writing very readable code. 

**Example**: let us assume that, internally we want the LFSR state stored as an integer with bits in reverse order to allow an efficient computation for the LFSR iterations, but, externally, i.e., when the state is read, we want to return the value with bits in correct order. We can define a private attribute `_state` for computations and a public attribute `state` that hides (i.e., invokes) a **getter method** that return the private attribute inverting the bit order.
```python
    @property
    def state(self):
        return reverse(self._state, len(self))
```
The same principle can be employed to let the attribute `state` invoke a **setter method** that perform the inverse operation:
```python
    @state.setter
    def state(self, state):
        self._state = reverse(state & self._statemask, len(self))
```
Note that the decorator now is `<attribute>.setter`.

Moreover, the setter method can be used to make attributes only readble and not writable. For example, we may want the LFSR length to be an only-read attribute.
```python
    @property
    def length(self):
        return self._length
    @length.setter
    def length(self, val):
        raise AttributeError('Denied')
```

Thanks to these feature we choose:
- `state` to be both redable and writable but setter and getter are defined in such a way to allow state to be internally stored with bits in reverse order.
- `poly` to be only redable and internally stored as integer.
- `length`, `output`
- `feedback` to be both redable and writable. The user can insert the state one bit at a time and we will see that some stream ciphers exploit this possibility.




```python
class LFSR(object):
    '''
    Class implementing a Linear Feedback Shift Register (LFSR)
    
    Attributes
    ----------
    poly: list of int,
        feedback polynomial expressed as list of integers corresponding to
        the degrees of the non-zero coefficients.
    state: int,
        state of the shift register.
    output: bool,
        last shift register output,
    length: int,
        length of the shift register as well as  maximum degree of the feedback
        polynomial.
    
    Methods
    -------
    run_steps(self, N=1)
        Execute multiple LFSR steps
    cycle(self, state=None)
        Execute a full cycle.
    '''
    
    def __init__(self, poly, state=None):
        '''
        Parameters
        ----------
        poly: list of int,
            feedback polynomial expressed as list of integers corresponding to
            the degrees of the non-zero coefficients.
        state: int, optional (default=None)
            shift register initial state.
            If None, state is set to all ones.
        '''
        length = max(poly)
        self._length = length
        self._poly = sparse_to_int(poly) >> 1 # p0 is omitted (always 1)
        
        self._statemask = (1 << length) - 1
        if state is None:
            state = self._statemask
        self.state = state
        
        self._outmask = 1 << (length-1)
        self._output = bool(self._state & self._outmask)
        
        self._feedback = parity(self._state & self._poly) 

    # ==== state ====
    @property
    def state(self):
        # state is re-reversed before being read
        return reverse(self._state, len(self))
    @state.setter
    def state(self, state):
        if not isinstance(state, int):
            raise TypeError('input type is not supported')
        # ensure seed is in the range [1, 2**m-1]
        if (state < 1) or (state > len(self)):
            state = 1 + state % (2**len(self)-2)
        self._state = reverse(state & self._statemask, len(self))

    # ==== length ====
    @property
    def length(self):
        return self._length
    @length.setter
    def length(self, val):
        raise AttributeError('Denied')
    
    # ==== poly ====
    @property
    def poly(self):
        return int_to_sparse((self._poly << 1) | 1)[::-1]
    @poly.setter
    def poly(self, poly):
        raise AttributeError('Denied')

    # ==== output ====
    @property
    def output(self):
        return self._output
    @output.setter
    def output(self, val):
        raise AttributeError('Denied')

    # ==== feedback ====
    @property
    def feedback(self):
        return self._feedback
    @feedback.setter
    def feedback(self, feedback):
        self._feedback = bool(feedback)
    
    
    def __str__(self):
        poly = ' + '.join([
            (f'x^{d}' if d > 1 else ('x' if d==1 else '1')) 
            for d in self.poly
        ])
        string = ', '.join([
            f'poly: "{poly}"',
            f'state: 0x{self.state:0{(self.length+1)//4}x}',
            f'output: {None if self.output is None else int(self.output)}'
        ])
        return string
        
    def __repr__(self):
        return f'LSFR({str(self)})'
    
    def __len__(self):
        return self.length
    
    def __iter__(self):
        return self
    
    def __next__(self):
        '''Execute a LFSR step and returns the output bit (bool)'''
        self._state = ((self._state << 1) | self._feedback) & self._statemask
        self._output = bool(self._state & self._outmask)
        self._feedback = parity(self._state & self._poly) 
        return self._output
    
    def run_steps(self, n=1):
        '''
        Execute multiple LFSR steps.
        
        Parameters
        ----------
        n: int, optional (default=1)
            number of steps to execute.
        
        Output
        ------
        list of bool (len=n),
            LFSR output bit stream.
        '''
        return [next(self) for _ in range(n)]
   
    def cycle(self, state=None):
        '''
        Execute a full LFSR cycle (LFSR.len steps).
        
        Parameters
        ----------
        state: int or list of int or bools, optional (default=None)
            shift register state. If None, state is kept as is.
        
        Output
        ------
        list of bool (len=2**myLFSR.len - 1),
            LFSR output bit stream.
        '''
        if state is not None:
            self.state = state
        return self.run_steps(n=int(2**len(self)) - 1)

```

Again, we check the implementation with the reference example. This time, we did not need to use a verbose flag to print the internal state, since all internal variable are accessible by the user.


```python
poly = [3, 1, 0]
state = 0x7
niter = 7

def print_lfsr(lfsr):
    print(f'{lfsr.state:0{len(lfsr)}b} ({lfsr.state:d}) ',
          f'{int(lfsr.output):d}  {int(lfsr.feedback):d}')

# create and instance of the LFSR
lfsr = LFSR(poly, state)

print('\n state   b fb')
print_lfsr(lfsr) # print initial state
for b in islice(lfsr, niter):
    print_lfsr(lfsr)
```

    
     state   b fb
    010 (2)  0  0
    001 (1)  1  1
    100 (4)  0  1
    110 (6)  0  1
    111 (7)  1  0
    011 (3)  1  1
    101 (5)  1  0
    010 (2)  0  0


Here is a check for the implementation of `__str__` and `__repr__` methods.
Remind that the goal of :
- `__str__` is to provide a readable string
- `__repr__` is to be unambiguous


```python
print(lfsr)   # print calls `__str__`
display(lfsr) # display calls `__repr__`
```

    poly: "x^3 + x + 1", state: 0x2, output: 0



    LSFR(poly: "x^3 + x + 1", state: 0x2, output: 0)


Here is a check for the implementation of the methods `cycle` and `run_steps`. Since the former calls the latter, we just check the former.


```python
cycle = lfsr.cycle()
print(cycle)
print(''.join(f'{int(b)}' for b in cycle))
```

    [True, False, False, True, True, True, False]
    1001110


## Berlekamp-Massey Algorithm

The Berlekamp–Massey algorithm is an algorithm that finds the shortest LFSR for a given binary sequence.
The input consists in the binary sequence and the output is the feedback polynomial characterizing the shortest LFSR able to generate that sequence.

**Pseudocode**:
> **Input** $b = [b_0, b_1, \dots, b_N]$ <br>
> $P(x) \leftarrow 1, \quad m \leftarrow 0 $ <br>
> $Q(x) \leftarrow 1, \quad r \leftarrow 1 $ <br>
> **for** $\tau = 0, 1, \dots, N-1 $ <br>
> $\qquad d \leftarrow \oplus_{j=0}^{m} p_j \otimes b[\tau-j] $ <br>
> $\qquad $ **if** $d = 1 $ **then** <br>
> $\qquad \qquad $ **if** $2m \le \tau $ **then** <br>
> $\qquad \qquad \qquad R(x) \leftarrow P(x) $ <br>
> $\qquad \qquad \qquad P(x) \leftarrow P(x) + Q(x)x^r $ <br>
> $\qquad \qquad \qquad Q(x) \leftarrow R(x) $ <br>
> $\qquad \qquad \qquad m \leftarrow \tau + 1 - m $ <br>
> $\qquad \qquad \qquad r \leftarrow 0 $ <br>
> $\qquad \qquad $ **else** <br>
> $\qquad \qquad \qquad P(x) \leftarrow P(x) + Q(x)x^r $ <br>
> $\qquad \qquad $ **endif** <br>
> $\qquad $ **endif** <br>
> $\qquad r \leftarrow r + 1 $ <br>
> **endfor** <br>
> **Output** $P(x) $ <br>

**Template**:
```python
def berlekamp_massey(b):
    ''' function docstring '''
    # algorithm implementation
    return poly
```

**Example**:

Input: $ b = [1, 0, 1, 0, 0, 1, 1, 1] $

| $\tau$ | $b_\tau$ | $d$ |      $P(x)$     | $m$ |    $Q(x)$   | $r$ |
|:------:|:--------:|:---:|:---------------:|:---:|:-----------:|:---:|
|   $-$  |    $-$   | $-$ | $ 1           $ | $0$ | $ 1       $ | $1$ |
|   $0$  |    $1$   | $1$ | $ 1 + x       $ | $1$ | $ 1       $ | $1$ |
|   $1$  |    $0$   | $1$ | $ 1           $ | $1$ | $ 1       $ | $2$ |
|   $2$  |    $1$   | $1$ | $ 1 + x^2     $ | $2$ | $ 1       $ | $1$ |
|   $3$  |    $0$   | $0$ | $ 1 + x^2     $ | $2$ | $ 1       $ | $2$ |
|   $4$  |    $0$   | $1$ | $ 1           $ | $3$ | $ 1 + x^2 $ | $1$ |
|   $5$  |    $1$   | $1$ | $ 1 + x + x^3 $ | $3$ | $ 1 + x^2 $ | $2$ |
|   $6$  |    $1$   | $0$ | $ 1 + x + x^3 $ | $3$ | $ 1 + x^2 $ | $3$ |
|   $7$  |    $1$   | $0$ | $ 1 + x + x^3 $ | $3$ | $ 1 + x^2 $ | $4$ |

Output: $ P(x) = 1 + x + x^3 $

We first define a functions to transform a bit sequence (list of bool, list of 1/0, or a string composed by only '1'/'0') into an integer.


```python
def bits_to_int(bits):
    ''' transform a bit sequence (str of 1/0 or list of bool) into an int '''
    integer = sum(1 << i for i, bit in enumerate(bits) if bool(int(bit)))
    return integer
```


```python
bits = '101011010111101'
integer = bits_to_int(bits)
bin(integer)
```




    '0b101111010110101'



We choose to implement the Berlekamp-Massy algorithm by storing polynoimals $P(x)$ and $Q(x)$ as integers. Instead, the bit sequence $b$ is kept as list/string and at each iteration is transformed to an integer.

Note that, from an implementation point of view, the main steps are:
- the computation of the discrepancy bit 
$$d \leftarrow \oplus_{j=0}^{m} p_j \otimes b[\tau-j]$$
This operation is very similar to the combination of state and polynominal coefficients in the computation of the LFSR feedback bit.
```python 
d = parity(P & bits_to_int(bits[t-m:t+1][::-1]))
```


- the update of the polynomial $P(x)$
$$P(x) \leftarrow P(x) + Q(x)x^r$$
```python 
P = P ^ (Q << r)
```
    - Assuming $Q(x) = \sum_{i=0}^l q_i x^i$, then $Q(x)x^r = \sum_{i=0}^l q_i x^{r+i}$. Therefore, if $Q(x)$ is represented as `[q0, q1, ..., ql]`, $Q(x)x^r$ is `[0, ..., 0, q0, q1, ..., ql]` where the number of 0s is $r$. With the integer representation, this is translated in a left shift of `Q` by `r` positions.
    - The $+$ operation in $GF(2^m)$ corresponds to the element-wise $+$ operation in $GF(2)$, that is the `xor`. With the integer representation the bit-wise xor is implemented with `^`.

As for the LFSR implemented as a generator, we make use of a `verbose` flag to print the internal variables for debugging purposes.


```python
def berlekamp_massey(bits, verbose=False):
    '''
    Find the shortest LFSR for a given binary sequence.
    
    Parameters
    ----------
    bits: list/string of 0/1,
        stream of bits.
        
    Return
    ------
    list of int,
        linear feedback polynomial expressed as list of the exponents related
        to the non-zero coefficients
    '''
    
    # variables initialization
    P, m = 0x1, 0
    Q, r = 0x1, 1
    
    if verbose:
        from math import log2
        lenP = 1 + int(log2(len(b)+1))
        print(f'  t b d {"poly":>{lenP}} m {"Q":>{lenP}} r')
        print(f'  - - - {P:{lenP}b} 0 {Q:{lenP}b} 1')
    
    for t in range(len(bits)):
        
        # compute discrepancy
        d = parity(P & bits_to_int(bits[t-m:t+1][::-1]))
        
        if d:
            if 2*m <= t: # A                
                P, Q = P^(Q<<r), P
                m, r = t+1-m, 0
            else: # B
                P = P^(Q<<r)
        r += 1
        
        if verbose:
            print(f'{t:3d} {bits[t]} {d:d} {P:{lenP}b} {m} {Q:{lenP}b} {r}')
            
    # convert poly from integer to list of degree of non-zero coefficients
    poly = int_to_sparse(P)[::-1]
    return poly
```

We can check the implementation by replicating the reference example.


```python
b = '1010011101'
# b = '11101001110100111010'
poly = berlekamp_massey(b, verbose=True)
print(poly)
```

      t b d poly m    Q r
      - - -    1 0    1 1
      0 1 1   11 1    1 1
      1 0 1    1 1    1 2
      2 1 1  101 2    1 1
      3 0 0  101 2    1 2
      4 0 1    1 3  101 1
      5 1 1 1011 3  101 2
      6 1 0 1011 3  101 3
      7 1 0 1011 3  101 4
      8 0 0 1011 3  101 5
      9 1 0 1011 3  101 6
    [3, 1, 0]


A more deep functioning assessment can consist in generating a bit sequence with an LFSR characterized by a primitive polynomial in the feedback, and check that the Berlekamp-Massey algorithm is able to infer the LFSR polyonimal.

Use LFSR with different length!


```python
# poly = [3, 1, 0]
# poly = [4, 1, 0]
# poly = [5, 2, 0]
# poly = [6, 1, 0]
# poly = [7, 1, 0]
# poly = [8, 4, 3, 2, 0]
# poly = [9, 4, 0]
# poly = [10, 3, 0]
poly = [11, 2, 0]
# poly = [16, 12, 3, 1, 0]

# create an instance of an LFSR
lfsr = LFSR(poly)
print(lfsr)

# generate a bit sequence with length of at least one cycle
bitstream = lfsr.cycle()
print(''.join([str(int(b)) for b in bitstream]))

# use the BM algorithm to infer the polynomial
bm_poly = berlekamp_massey(bitstream)
print(bm_poly)
```

    poly: "x^11 + x^2 + 1", state: 0x002, output: 0
    1000000000010101010101110111011100100111001011110010111001111011000011000111111000001100111111100001100110101001011011101000111001000001101000000011101010101100010001001010000010111010101110010001001111010111100111011001011000111101101011001001000111101000001100010101001010001000010000010101111111011100110010011110000101100101011100001000110101000001110111111100100110011110100101100100001001011111101000110011101011010011101101111001001001101000010110111111011011001101101101001001001000010111101010001100010000011111010101100111011100001101100000011100000000110101010100100010001011111010111001101110010110110000100100101011110100010011011111010010011000101111000001001101010111100010001100000101001010100010111011111011000110010010100101111011110110010011011010000111000101011000001000111111101011001100010010110101111011011100111000110100111110001011001111101101001100011101001010011101000101100010100011111011111001101100111100011110011111001111001101001101001000011101000000110111111110001100110000011110000000110000000001111111111100110011001011010010111000101110010100011010001010010000010000101010000001000101010111110111011001110010010110000101110000001001111111101001100110111100001110011111100101100110100011110001010011010111010010001101111010110110010001110000101001111110111100110001100001111100000011001010101101000100011101011111001000110010111110000100110000001011010101000111011101011000110111000001110010101001111011101001101100010110110101110001110110000011011010101101101110110110110001110001111100101001100001000011111101010011001000101101000001001000000010111111111011001100111000011110010101100101110110100011011011111000111001101011000011101101010011100010000110101111110001001100101000011110111111001110011000011010010101101111011100011000110101101011011100010010011111010000110010000001111010101001101110111100011011001010010010111010000100111010100001101110101001001110111101001110011101001111001000011000010101101010001001000101000010100010101110101110110001001110000010110000000100101010100001000100000010100000000
    [11, 2, 0]


## Berlekamp-Massey Streaming Algorithm

An algorithm is streaming when the input is a stream of symbols. That means the input is a sequence of items and the algorithm is applied to only one item (or a group of successive items) at a time.

Since Berlekamp-Massey Algorithm may be applied to very long sequence of bits,  it may be convenient to implement it as a streaming algorithm. The streaming implementation allows for lower memory requirements as you need to store only the bits of the sequence that are significant for the identification of the LFSR polynomial.

The streaming algorithm must:
- take one bit at a time as input
- for each input bit, it must return the polynomial corresponding to the shortest LFSR that can generate the bit sequence observed up to the last input bit.

Since streaming algorithms needs an internal state to be updated at each new incoming symbol of the stream, they cannot be implemented a normal functions. We can implement it as a class. 

**Template**:
```python
class BerlekampMassey():
    
    def __init__(self):
        # do stuff
        self.poly = ...
        
    def __call__(self, bit):
        # do stuff
        return self.poly
```
The special method `__call__` makes the class callable, that means you can call
it as a function.
```python
berlekamp_massey = BerlekampMassey()

for bit in bit_stream:
    poly = berlekamp_massey(bit)
```


Below is the implementation of the Berlekamp-Massey algorithm as a streaming algorithm.

Unlike the batch implementation, here the bits come one at a time and therefore the algorithm needs to internally store the bit sequence by appending each new incoming bit. We choose to store the bit sequence `bits` as an integer, and the new incoming bit `bit` is appended as least significant bit (LSB):
```python
self.bits = (self.bits << 1) | int(bit)
```

Moreover, one can notice that not all bits of the sequence are employed to determine the polynomial so we can drop them. It can be proven that the drop operation can be performed when $Q(x)$ is updated and $r$ is reset, by selection only the last $m$ bits of the sequence.
```python 
self.bits &= (1 << self.m) - 1
```


```python
class BerlekampMassey():
    '''
    Berlekamp-Massey Algorithm. 
    The algorithm finds the shortest LFSR for a given binary sequence.
    
    This class returns a function whose call method takes one bit at a time as 
    input and returns the feedback polynomial of the shortest LFSR capable of 
    generating the sequence of all bits received up to the last one.
    
    Attributes
    ----------
    poly: list of int,
        feedback polynomial expressed as list of integers corresponding to
        the degrees of the non-zero coefficients.
    
    Methods
    -------
    __call__(bit):
        update and return the polynomial of the shortest LFSR for the observed 
        bit sequence.
    discrepancy():
        compute the discrepancy bit that is 0 (False) if the polynomial `poly` 
        explain the observed bit sequence, or 1 (True) otherwise.
    '''
    
    def __init__(self):
        ''' class constructor '''
        # init algorithm's internal variables
        self.P, self.m = 0x1, 0
        self.Q, self.r = 0x1, 1
        # init an empty sequence of bits
        self.bits = 0 # self.bits = []
        self.tau = 0
    
    # ==== poly ====
    @property
    def poly(self):
        return int_to_sparse(self.P)[::-1]
    @poly.setter
    def poly(self, val):
        raise AttributeError('Denied')

    def discrepancy(self):
        b = self.bits & ((1 << (self.m + 1)) - 1)
        d = parity(self.P & b)
        return d
    
    def __call__(self, bit):
        '''
        Update the feedback polynomial characterizing the shortest LFSR capable 
        of generating the sequence of all bits received.
    
        Parameters
        ----------
        bit: bool, int, or 0/1 str
            input bit
        
        Return
        ------
        poly: list of int,
            feedback polynomial expressed as list of integers corresponding to
            the degrees of the non-zero coefficients.
        '''
        # append
        self.bits = (self.bits << 1) | int(bit)  # self.bits.append(bit)
        if self.discrepancy():
            if 2*self.m <= self.tau: # A                
                self.P, self.Q = self.P ^ (self.Q << self.r), self.P
                self.m, self.r = self.tau + 1 - self.m, 0
                self.bits &= (1<<self.m) - 1  # self.bits = self.bits[-self.m:]
            else: # B
                self.P = self.P ^ (self.Q << self.r)
        self.r += 1
        self.tau += 1
        
        return self.poly
```

Again, we can assess the implementation of the algorithm by generating a sequence with an LFSR featuring a primitive polynomial and check whether the algorithm is able to determine the polynomial.


```python
# poly = [3, 1, 0]
poly = [11, 2, 0]

lfsr = LFSR(poly)

berlekamp_massey = BerlekampMassey()
for bit in islice(lfsr, 2**len(lfsr)-1):
    poly = berlekamp_massey(bit)
poly
```




    [11, 2, 0]



A deeper check may be performed by trying to determine the polynomial pf the LFSR generating the sequence but trying with any possible initial state.


```python
import tqdm # very useful module to implement progress bars
```


```python
# poly = [3, 1, 0]
# poly = [9, 4, 0]
poly = [11, 2, 0]

# iter over every possible initial state for the LFSR, i.e., from 1 to 2^m-1
for state in tqdm.tqdm(range(1, 2**max(poly)-1), ncols=79):
    
    # create and instance of the LFSR with the initial state 
    lfsr = LFSR(poly, state)
    
    # create an instance of the BM algorithm
    berlekamp_massey = BerlekampMassey()
    
    # iter over the bit sequence and apply the algorithm
    for bit in islice(lfsr, 2**max(poly)):
        poly_ = berlekamp_massey(bit)
        
    # check 
    if poly_ != poly:
        print(f'algorithm failed for initial state 0x{state:X}.')
        break
        
else:
    print('algorithm succeded for all possible initial states.')
```

    100%|██████████████████████████████████████| 2046/2046 [01:30<00:00, 22.59it/s]

    algorithm succeded for all possible initial states.


    


78fb81c3820d3dd9b0d358ba03531d563efdc23e
ipynb
Jupyter Notebook
3.lfsr/LFSR.ipynb
marchioa/data-security
2bc594c85bfbf0636df9d5a5f688e9c8b9fe12e3
MIT
3.lfsr/LFSR.ipynb
marchioa/data-security
2bc594c85bfbf0636df9d5a5f688e9c8b9fe12e3
MIT
3.lfsr/LFSR.ipynb
marchioa/data-security
2bc594c85bfbf0636df9d5a5f688e9c8b9fe12e3
MIT
Qwen/Qwen-72B
1. YES
2. YES