The challenge gives us a Fibonacci LFSR, with unknown taps, but we are provided 2048 bits of the output of LFSR. We have to first recover the taps, then from the taps we will recover the initial value.

To recover the taps from the given output, we can use the Berlekamp Massey algorithm. We are provided with 2048 bits, which is more than enough (we need more than $384 * 2 = 768$ bits). Sage has its own implementation of this algorithm, more information can be found here.

Instead of this approach, as we are provided with $2n$ bits of output ($n = 384$), we can take on this challenge using an approach using an approach from an angle of linear equations. This is inspired from this blog entry (look at the “stupid” approach) and solving system of linear equations using matrices. Kudos to yassinebelarbi for helping me, I was stuck on this for 2 weeks.

First, xor operation is addition in GF(2). We can think of the Fibonacci LFSR, which uses the xor operation, as a linear combination of the bits. Denote the coefficients (taps) as $c_0, c_1, c_2, …, c_{n-1}$. Denote the state as $a_0, a_1, a_2, …, a_{n-1}$.

Suppose the tap values are $(0, 5, 15, 16)$, $c_0 = 1, c_5 = 1, c_{15} = 1, c_{16} = 1$, and the other coefficient values are $0$.

The output of the LFSR is the linear combination $c_0 a_0 + c_1 a_1 + c_2 a_2 + … + c_{n-1} a_{n-1}$ in $GF(2)$. The output of the linear combination will be $a_n$. Denote $S_k = (a_k, a_{k + 1}, a_{k + 2}, …, a_{k + n})$, the result of the linear combination of $S_k$ and the coefficients is $a_{k + n + 1}$.

Since we know $2n$ bits, we know the output of the LFSR for the first $n$ iterations. More specifically, we know the for each state $S_k$ $(0 < k < n)$, the corresponding output $a_{k + n}$. Refer to this link for the following matrix construction.

We can construct a matrix $A = (S_0, S_1, …, S_{n - 1})$ (the Ss are the rows of the matrix). The resultant vector is $a_n, a_{n+1}, …, a_{2n -1}$. Again, refer to the link for how to solve this problem - it is simple linear algebra mathematics. We should be able to recover the coefficients (taps). Refer to this page about Companion Matrix for more information on the “meaning” of the coefficients. TLDR, $c_0$ corresponds to the coefficient of the lowest degree component, and $c_{n-1}$ corresponds to the coefficient of the highest degree component.

Sage implementation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

output = "00101101100100011101101110000000101100001110110001011110110011011011011011001111111101111101011111010000001011000011111010110011010111000100110010101100101000011000111010110011110011001001000101010010011100001011100101010000111100001010101001100010101110110110111000111010000011000000000010110101011011011111100100010100101100110000111010001101010111010010011110000110110010100010011101111111101001111101001100110110111110100101000010000000000010000000100100001110001110100001100111111010110010010110000000001110101101011001011011001101001010000000100010011001000100111000110100101011101001100110010110010011100010011101010000010101010001111110101010111001111000111111011011101011100011010101011010110001101110110101111001101110010100010110100100000000100000111000110010101100111111101001100111110000101011100111110010001111100000110110110110010101101110111100101010011000100110001000110011111110011110111001100110001011111111111011010010010110010011101011001010010101011001010100110001001111111000110011010101011011010110000101011011001100100000000000011111110101001100011001110000101110001101011011110000100111100010111110101101101010100110100010011001000110010100000001110010011001111110010001110110110100000001110110111011010100100111011011100001110100001110101111010100110101110100111111101100011010110011111111100011101110001111011100000000100110010111100100110100011000110110001001011111110000111001110010010111101000100011100101001111101011011011111011100011000001011011101111101111011011101011101110011011111001011101101100111111000000100000100011111101110100110100000111010001101100011000100000100010010111101111011000001011110101001010001000001001000111101110101001010110101101101010001110000001000001011110001011100011011010010110011011111111110010111101011111111101111110001100101001011010110001000001000111000101010100001000000101000100011100111001100111010010100101100100101001010100001100111001011111110011000000110000110101101000001110100101110011110001111000001111000100111000111011010100101000001101111100011011000111010001110110"

# Convert the given output to a stream
stream = [int(x) for x in [*output]]

# Construct matrix to obtain the tap values
mat = []
for i in range(384):
    mat.append(stream[i:i+384])

mat = Matrix(GF(2), mat)

# The vector is the output of each row in the matrix, 
# which is the next 384 bits in the given output
vec = vector(GF(2), stream[384 : 384 * 2])

print(mat.inverse() * vec)

With the taps figured out, now we can use the Matrix form of Fibonacci LFSR. $k = 16 \times 48$, as before the given output, there is $16 \times 48$ bits produced in the init method.

The following Sage implementation uses the Berlekamp Massey algorithm, and should detail how to get the initial bit string (flag) with the given tap values.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48


from Crypto.Util.number import long_to_bytes
from sage.matrix.berlekamp_massey import berlekamp_massey

# Output from the given output 
output = "00101101100100011101101110000000101100001110110001011110110011011011011011001111111101111101011111010000001011000011111010110011010111000100110010101100101000011000111010110011110011001001000101010010011100001011100101010000111100001010101001100010101110110110111000111010000011000000000010110101011011011111100100010100101100110000111010001101010111010010011110000110110010100010011101111111101001111101001100110110111110100101000010000000000010000000100100001110001110100001100111111010110010010110000000001110101101011001011011001101001010000000100010011001000100111000110100101011101001100110010110010011100010011101010000010101010001111110101010111001111000111111011011101011100011010101011010110001101110110101111001101110010100010110100100000000100000111000110010101100111111101001100111110000101011100111110010001111100000110110110110010101101110111100101010011000100110001000110011111110011110111001100110001011111111111011010010010110010011101011001010010101011001010100110001001111111000110011010101011011010110000101011011001100100000000000011111110101001100011001110000101110001101011011110000100111100010111110101101101010100110100010011001000110010100000001110010011001111110010001110110110100000001110110111011010100100111011011100001110100001110101111010100110101110100111111101100011010110011111111100011101110001111011100000000100110010111100100110100011000110110001001011111110000111001110010010111101000100011100101001111101011011011111011100011000001011011101111101111011011101011101110011011111001011101101100111111000000100000100011111101110100110100000111010001101100011000100000100010010111101111011000001011110101001010001000001001000111101110101001010110101101101010001110000001000001011110001011100011011010010110011011111111110010111101011111111101111110001100101001011010110001000001000111000101010100001000000101000100011100111001100111010010100101100100101001010100001100111001011111110011000000110000110101101000001110100101110011110001111000001111000100111000111011010100101000001101111100011011000111010001110110"

# Convert the given output to a stream
stream = [int(x) for x in [*output]]
# print(berlekamp_massey(stream))

# Taps recovered from Berlekamp Massey algorithn
taps = [16, 15, 6, 0]

# Coefficients for the last row in the matrix representation of the Fibonacci LFSR
# See the matrix at https://en.wikipedia.org/wiki/Linear-feedback_shift_register#Matrix_forms
coeffs = [0 for i in range(384)]

for tap in taps: 
    coeffs[tap] = 1

# Constructing the Fibonacci LFSR
mat = []

for i in range(383):
    temp = [0 for i in range(384)]
    temp[i + 1] = 1 
    mat.append(temp) 

mat.append(coeffs)
mat = Matrix(GF(2), mat)

# The vector is the first 384 bits from the output
# Again, see the matrix form of LFSR in the Wikipedia link above, 
# this vector is the LHS vector
vec = vector(GF(2), stream[:384])

# k = 16 * 48 (from LFSR matrix form), and the init method of LFSR
mat = mat ^ (16 * 48)

# Obtain the flag, this is similar to solving linear equations using matrix
# See this https://www.mathsisfun.com/algebra/systems-linear-equations-matrices.html, 
# We are solving AX = B, and X = A ^ -1 * B
flag = list(mat.inverse() * vec)
flag = [int(x) for x in flag]
flag = ''.join(str(x) for x in flag)

print(long_to_bytes(int(flag, 2)))

One thing to note, and this burns A LOT of my time in solving, is that during testing, we cannot just pick any arbitrary value of taps and expect the above to work (even with the BM algorithm). The taps must be chosen such that the polynomial is minimal. The definition of a minimal polynomial is a polynomial of an algebraic number is the unique irreducible monic polynomial of smallest degree with 1 as the leading coefficient.

The method ONLY works on the given output, and I was too dumb to realize this.