Predicting Fibonacci LFSRs over Integer Residue Rings Using Multi-Segment Truncated Sequences

'''
This code is used to recover the truncated ZUC's driving sequence by multi-segment truncated sequences.

Suppose that the modulus, the coefficients and the initial state are unknown, given the first N truncated digits which consist of the middle segment(s) of the sequence, we need to recover the modulus, the coefficients and the initial state.

'''
import time

# Search annihilating polynomials over Z/(m) of the sequence
def Search_linear_relations(truncated_a, v, r, t):
    l = len(v)
    L = Matrix(ZZ, r + t, r + t)
    L[0 : t, 0 : t] = identity_matrix(t) * 2 ** (v[l - 1] - v[1])
    L[t : r + t, t : r + t] = identity_matrix(r) * 2 ** (v[l - 2] - v[1])
    for i in range(r):
        for j in range(t):
            for w in range((l - 1) // 2):
                L[i + t, j] += 2 ** (v[2 * w + 1] - v[1]) * truncated_a[2 * w + 1, i + j]
    reduced_L = L.BKZ(algorithm = 'NTL', fp = 'fp', block_size = 20)
    C = reduced_L[0 : r, t : r + t] / 2 ** (v[l - 2] - v[1])

    total_eta_m = 0
    for j in range(r):
        s_m = 0
        if C[j] != 0:
            for i in range(t):
                if (C[j] * a[i : i + r]) % m == 0:
                    s_m += 1
            if s_m == t:
                total_eta_m += 1
    print("Number of annihilating polynomials = ", total_eta_m)
    return C

# Determine whether c is the n-th power of some k-bit prime
def is_nth_power(c, n, k):
    if ZZ(c).is_perfect_power():
        base, exp = ZZ(c).perfect_power()
        if base.is_prime() and base.nbits() == k and exp == n:
            return base
        else:
            return False
    else:
        return False

# Recover the modulus
def recover_the_modulus(C, n, k):
    R = PolynomialRing(ZZ, 'x')
    (x,) = R._first_ngens(1)
    r = C.ncols()
    f = []
    fi = C[0, 0] * x**0
    for i in range(1, r):
        fi += C[0, i] * x**i
    f.append(fi)
    fi = C[1, 0] * x**0
    for i in range(1, r):
        fi += C[1, i] * x**i
    f.append(fi)
    m_power = f[0].resultant(f[1])

    flag = 2
    recovered_m = is_nth_power(m_power, n, k)

    while recovered_m == False and flag < r:
        fi = C[flag, 0] * x**0
        for i in range(1, r):
            fi += C[flag, i] * x**i
        f.append(fi)
        for i in range(flag):
            res = f[flag].resultant(f[i])
            m_power = gcd(m_power, res)
        flag += 1
        recovered_m = is_nth_power(m_power, n, k)

    if recovered_m == False:
        Flag = 0
        print("Fail to recover m!")
        #print("The wrong m == ", recovered_m, "\n")
    else:
        Flag = 1
        print("Success to recover m by ", flag, " polynomials!")
        print("m = ", recovered_m, "\n")

    return recovered_m

# Recover the coefficients
def recover_the_coefficients(C, m, n):
    R = PolynomialRing(Integers(m), 'x')
    (x,) = R._first_ngens(1)
    r = C.ncols()

    for i in range(r):
        mark = r - 1
        while C[i, mark] == 0:
            mark -= 1
        for j in range(r):
            C[i, j] = (C[i, j] * C[i, mark]**(-1)) % m

    f = []
    fi = ZZ(C[0, 0]) * x**0
    for i in range(1, r):
        fi += ZZ(C[0, i]) * x**i
    f.append(fi)
    fi = ZZ(C[1, 0]) * x**0
    for i in range(1, r):
        fi += ZZ(C[1, i]) * x**i
    f.append(fi)
    g = f[0].gcd(f[1])

    flag = 2
    while g.degree() > n and flag < r:
        fi = ZZ(C[flag, 0]) * x**0
        for i in range(1, r):
            fi += ZZ(C[flag, i]) * x**i
        f.append(fi)
        for i in range(flag):
            g = g.gcd(f[flag])
        flag += 1

    if g.degree() != n:
        Flag = 0
        print("Fail to recover f!")
        #print("The wrong f = ", g, "\n")
    else:
        Flag = 1
        print("Success to recover f by ", flag, " polynomials!")
        print("f = ", g, "\n")

    return g

def recover_the_initial_state(m, n, f, truncated_a, k, v, d):
    l = len(k)
    dim = d * (l + 1) // 2 + 1
    recovered_a = vector(ZZ, n)
    Q = Matrix(ZZ, n, n)
    Q_power = identity_matrix(n)
    L = Matrix(ZZ, dim, dim)

    Q[0, n - 1] = (-f.list()[0]) % m
    for i in range(1, n):
        Q[i, i - 1] = 1
        Q[i, n - 1] = (-f.list()[i]) % m

    for i in range(1, n):
        Q_power = (Q_power * Q) % m
    for i in range(n, d):
        Q_power = (Q_power * Q) % m
        b = 0
        L[1 : n + 1, i + 1] = Q_power[0 : n, 0]
        for j in range(int((l - 1) / 2)):
            L[j * n + d + 1 : j * n + d + n + 1, i + 1] = Q_power[0 : n, 0] * 2**v[2 * (j + 1)]
            b += 2 ** v[2 * j + 1] * truncated_a[2 * j + 1, i]
            for w in range(n):
                b -= Q_power[w, 0] * 2 ** v[2 * j + 1] * truncated_a[2 * j + 1, w]
            b = b % m
        L[0, i + 1] = b + 2**(k[0] - 1)

    for i in range(n + 1):
        L[0, i] = 2**(k[0] - 1)
    for i in range(d + 1, dim):
        L[0, i] = 2**(k[0] - 1)
    L[1 : n + 1, 1 : n + 1] = identity_matrix(n)
    L[n + 1: d + 1, n + 1: d + 1] = identity_matrix(d - n) * m
    for i in range(int((l - 1) / 2)):
        L[i * n + d + 1: i * n + d + n + 1, i * n + d + 1: i * n + d + n + 1] = identity_matrix(n) * 2**(k[0] - k[2 * (i + 1)])
        L[i * (d - n) + d + (l - 1) // 2 * n + 1: i * (d - n) + d + (l - 1) // 2 * n + d - n + 1, n + 1 : d + 1] = identity_matrix(d - n) * (-2**v[2 * (i + 1)])
        L[i * (d - n) + d + (l - 1) // 2 * n + 1: i * (d - n) + d + (l - 1) // 2 * n + d - n + 1, i * (d - n) + d + (l - 1) // 2 * n + 1: i * (d - n) + d + (l - 1) // 2 * n + d - n + 1] = identity_matrix(d - n) * 2**(k[0] - k[2 * (i + 1)])

    reduced_L = L.BKZ(algorithm = 'NTL', fp = 'fp', block_size = 20)

    flag = 0
    while reduced_L[flag, 0].abs() != 2**(k[0] - 1):
        flag += 1
    if reduced_L[flag, 0] == -2**(k[0] - 1):
        for i in range(0, n):
            recovered_a[i] = reduced_L[flag, i + 1] + 2**(k[0] - 1)
            for j in range(int((l - 1) / 2)):
                recovered_a[i] += truncated_a[2 * j + 1, i] * 2**v[2 * j + 1]
                recovered_a[i] += (reduced_L[flag, j * n + d + 1 + i] + 2**(k[0] - 1)) // 2**(k[0] - k[2 * (j + 1)]) * 2**v[2 * (j + 1)]
            recovered_a[i] = recovered_a[i] % m
    if reduced_L[flag, 0] == 2**(k[0] - 1):
        for i in range(0, n):
            recovered_a[i] = -reduced_L[flag, i + 1] + 2**(k[0] - 1)
            for j in range(int((l - 1) / 2)):
                recovered_a[i] += truncated_a[2 * j + 1, i] * 2**v[2 * j + 1]
                recovered_a[i] += (-reduced_L[flag, j * n + d + 1 + i] + 2**(k[0] - 1)) // 2**(k[0] - k[2 * (j + 1)]) * 2**v[2 * (j + 1)]
            recovered_a[i] = recovered_a[i] % m

    return recovered_a


n = 16
m = 2147483647
# Exp.1
r = 70
# Exp.2
#r = 105
# Exp.3
#r = 104
# Exp.4
#r = 90
t = r
N = r + t - 1
a = vector(ZZ, N)

a[0] = 5113033
a[1] = 101513454
a[2] = 342517303
a[3] = 602278014
a[4] = 545182365
a[5] = 2109832480
a[6] = 83899935
a[7] = 1353806595
a[8] = 56274808
a[9] = 220365776
a[10] = 1900047307
a[11] = 1287703311
a[12] = 754976298
a[13] = 344592060
a[14] = 2072123851
a[15] = 1906585150
for i in range(n, N):
    a[i] = (2^15 * a[i - 1] + 2^17 * a[i - 3] + 2^21 * a[i - 6] + 2^20 * a[i - 12] + 257 * a[i - 16]) % m

# Exp.1 and Exp.2
l = 3
# Exp.3 and Exp.4
#l = 5
# Exp.1
k = vector([9, 17, 5])   # l = 3, r = 70, d = 30
# Exp.2
#k = vector([12, 14, 5])   # l = 3, r = 105, d = 36
# Exp.3
#k = vector([6, 3, 3, 14, 5])   # l = 5, r = 104, d = 30
# Exp.4
#k = vector([6, 2, 3, 15, 5])   # l = 5, r = 90, d = 31

v = vector(ZZ, l)
for i in range(l):
    for j in range(i):
        v[i] += k[j]

truncated_a = Matrix(ZZ, l, N)
for i in range(l):
    for j in range(N):
        truncated_a[i, j] = (a[j] >> v[i]) % 2 ** k[i]

print("Number of segments l = ", l)
print("Bit-length of each segment k = ", k, "\n")
print("******************************************************")
print("Step 1: Search annihilating polynomials over Z/(m) of the sequence")
Tstart = time.perf_counter()
C = Search_linear_relations(truncated_a, v, r, t)
Tend = time.perf_counter()
elapsed = Tend - Tstart
print("Wall Time: ",elapsed ," s\n")

print("******************************************************")
print("Step 2: Recover the modulus")
Tstart = time.perf_counter()
recovered_m = recover_the_modulus(C, n, 31)
Tend = time.perf_counter()
elapsed = Tend - Tstart
print("Wall Time: ",elapsed ," s\n")

print("******************************************************")
print("Step 3: Recover the coefficients")
Tstart = time.perf_counter()
recovered_f = recover_the_coefficients(C, recovered_m, n)
Tend = time.perf_counter()
elapsed = Tend - Tstart
print("Wall Time: ",elapsed ," s\n")
'''
R = PolynomialRing(Integers(m), 'x')
(x,) = R._first_ngens(1)
f = x**16 - (2**15 * x**15 + 2**17 * x**13 + 2**21 * x**10 + 2**20 * x**4 + 2**8 + 1)
'''

# Exp.1
d = 30
# Exp.2
#d = 36
# Exp.3
#d = 30
# Exp.4
#d = 31

print("******************************************************")
print("Step 4: Recover the initial state")
Tstart = time.perf_counter()
recovered_a = recover_the_initial_state(recovered_m, n, recovered_f, truncated_a, k, v, d)
print("a = ", recovered_a)
Tend = time.perf_counter()
elapsed = Tend - Tstart
print("Wall Time: ",elapsed ," s\n")